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This  dissertation  investigates  various  approaches  to  the  processing  of  narrowband 
frequency  data  from  large  aperture  acoustic  vertical  line  arrays  with  on  the  order  of 
one  hundred  equally  spaced  sensors.  When  little  is  known  about  the  oceanic 
environment,  beamforming  is  performed  to  detect  signals  and  estimate  the  vertical 
arrival  structure  of  the  pressure  field.  The  attention  is  focussed  on  using  the  adaptive 
Minimum  Variance  processor.  The  issues  of  signal  cancellation  due  to  correlated 
arrivals  and  mismatch  due  to  wavefront  curvature  are  examined.  In  practice,  planar 
wavefront  replica  vectors  do  not  result  in  significant  mismatch  for  a  900  m  aperture. 
Processing  subapertures  is  proposed  because  of  the  adaptive  processor  requirement  of 
having  numerous  data  snapshots  to  estimate  and  invert  the  covariance  matrix,  while 
the  pressure  field  may  not  be  stationary.  Full  aperture  and  subaperture  processing 
techniques  are  validated  with  simulated  data  created  by  a  normal  mode  model,  and 
then  used  on  200  Hz  data  collected  in  the  North-East  Pacific  by  the  Marine  Physical 
Laboratory  in  September  1987.  Multipath  arrivals  are  found  to  illuminate  only  parts  of 
the  120  sensor  array,  cut  for  100  Hz  and  900  m  long,  thus  indicating  that  the  wavefield 
can  be  highly  inhomogeneous  with  depth..  The  processing  gain  obtained  with  the 
adaptive  processor  is  found  lower  than  with  the  conventional  beamformer  because  of 
the  reduction  of  the  aperture  after  spatial  smoothing,  residual  signal  cancellation  and 
mismatch,  but  this  is  compensated  for  by  a  much  larger  resolution  for  high  level 
arrivals.  Beamforming  simulated  or  real  data  produces  a  gain  6  to  7  dB  lower  than  the 
array  gain  expected  with  matched  field  processing.  Matched  field  processing  on  the 
real  data  is  very  succesful  because  of  their  good  quality,  adequate  knowledge  of  the 
environment,  and  proper  acoustic  modeling.  The  200  Hz  source  is  correctly  localized, 
165  km  from  the  array.  Array  gain  estimates  indicate  that  matched  field  processing  is  a 
better  beamforming  approach  for  signal  detection  than  is  conventional  beamforming. 


Introduction 

Vertical  line  arrays  with  equally  spaced  sensors  are  systems  commonly  used  by 
the  scientific  community  as  experimental  tools  to  measure  ambient  noise  as  well  as 
long  range  propagation  signals  [Anderson,  1974,  Anderson,  1979,  Kewley,  1984, 
Dosso,  1987,  Sen,  1988,  Sotirin,  1989].  At  first  with  moderate  lengths  on  the  order 
of  a  few  hundred  meters,  array  apertures  recently  have  been  increased  to  provide 
better  coverage  of  the  water  column  and  to  enable  work  at  lower  frequencies  (from 
tens  to  hundreds  of  Hertz).  An  example  of  a  large  and  well  filled  acoutic  array  is  the 
Marine  Physical  Laboratory  (MPL)  digital  array,  a  120  sensor  aperture  of  900  m  cut 
for  100  Hz  [Sotirin,  1989].  This  dissertation  investigates  various  ways  to  process 
narrowband  signal  data  from  large  vertical  line  arrays  in  a  realistic  oceanic 
environment,  from  beamforming  which  theoretically  requires  little  or  only  very 
limited  knowledge  of  the  environment  to  matched  field  processing  where  the 
environment  has  to  be  fully  characterized. 

The  data  from  an  array  such  as  the  MPL  array  are  processed  using  spectral 
estimation  techniques  to  obtain  the  pressure  field  angular  spectrum  and  a  picture 
of  directionality.  Peaks  in  the  angular  spectrum  allow  the  identification  of  the 
direction  of  arrival  of  underwater  signals  with  enhanced  signal-to-noise  ratio.  The 
conventional  beamformer  is  commonly  used  due  to  its  computationally  efficient 
implementation  with  Fast  Fourier  Transforms  (FFTs)  and  its  robustness.  Its 
resolution  is  characterized  by  the  ratio  of  the  acoustic  wavelength  to  the  overall 
aperture  length  so  that  improvement  in  resolution  at  low  frequency  (or  large 
wavelength)  can  only  be  realized  with  large  arrays.  Since  feasibility  and  cost 
considerations  generally  put  a  limit  on  the  spatial  coverage,  processing  techniques 
with  increased  resolution  capabilities  are  very  attractive.  Adaptive  methods,  such  as 
the  Minimum  Variance  processing  technique  [Capon,  1969],  are  more  sophisticated 
and  make  specific  assumptions  on  the  wavefield  to  produce  such  resolution 
improvements  at  high  signal-to-noise  ratio.  A  benefit  in  using  adative  techniques  is 
the  ability  to  reject  strong  interferences  and  allow  the  detection  of  weaker  signals. 
An  increased  resolution  allows  the  separation  of  the  multipath  arrivals  impinging 
on  the  vertical  line  array,  thus  making  possible  a  detailed  study  of  the  complex 
propagation  physics.  Since  the  model  assumptions  made  by  the  adaptive  processor 
often  do  not  hold  in  actual  operation  [Vural,  1979],  the  influence  of  the  complex 
oceanic  environment  on  the  processing  is  addressed  here.  The  pressure  field  may 
be  strongly  influenced  by  the  sound  speed  profile  and  the  ocean  boundaries,  since 
large  aperture  arrays  have  lengths  on  the  order  of  the  water  depth.  The  variations 
of  sound  speed  with  depth  result  in  waveguide  propagation  where  the  arrivals  have 
curved  wavefronts.  The  bottom  and  surface  boundaries  combined  with  the 
waveguide  structure,  produce  multipath  arrivals  or  spatially  correlated  arrivals. 
Another  problem  with  the  adaptive  methods  is  their  requirement  for  large  data 
sample  size  to  estimate  the  array  covariance  matrix  and  obtain  a  statistically  stable 
estimate  of  the  power  spectrum.  This  can  be  troublesome  if  the  field  is  not 
stationary  such  as  in  the  case  of  a  moving  source. 

Chapter  1  discusses  the  array  processing  hypothesis,  and  presents  the 
conventional  and  adaptive  Minimum  Variance  processors.  The  general  beamforming 
procedure  is  decomposed  into  two  components,  a  processing  structure  and  a  signal 
model  which  assumes  uncorrelated  arrivals  with  given  wavefronts  (generally 
planar).  Implementation  issues  for  the  Minimum  Variance  processor  are  discussed. 
Chapter  2  investigates  two  aspects  of  sound  propagation  in  the  ocean:  the  problem 
of  multipaths  and  wavefront  curvature.  The  impact  of  these  environmental  factors 
on  the  conventional  beamformer  and  the  Minimum  Variance  processor  are  studied 
and  evaluated,  and  correcting  actions  examined.  Spatial  smoothing  techniques 
[Shan,  1985]  are  discussed  and  curved  wavefront  replica  vectors,  derived  with  a  ray 
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geometric  approach,  are  introduced.  Using  these  preliminary  results,  two  processing 
strategies  with  the  Minimum  Variance  processor  are  proposed  in  Chapter  3  and 
studied  through  simulations.  Full  array  aperture  and  array  subaperture  processing 
is  envisioned  with  curved  and  planar  wavefronts.  The  choice  of  the  most  adequate 
scheme  is  based  on  beamforming  a  3000  m  or  a  900  m  aperture  where  pressure 
field  data  at  a  particular  frequency  are  generated  by  an  acoustic  normal  mode 
model.  After  this  validation  on  synthetic  data,  these  processing  techniques  are  used 
in  Chapter  4  on  real  data,  collected  at  sea  by  the  MPL  digital  array  during  the 
September  1987  Vertical  Line  Array  (VLA)  Experiment  in  the  North-East  Pacific.  200 
Hz  signals  are  processed  and  their  arrival  structure  interpreted  with  ray  and  normal 
mode  acoustic  models.  Performance  aspects  (processing  gain  and  resolution),  are 
considered  for  various  input  signal-to-noise  ratios  in  the  case  of  the  synthetic  and 
the  real  data.  Finally,  Chapter  5  presents  matched  field  processing  which  is  a 
generalized  form  of  beamforming  where  a  global  knowledge  of  the  environment  is 
required.  Matched  field  processing  is  performed  on  the  same  data  set  as  in  Chapter 
4.  Localization  and  processing  gain  performance  are  examined  and  sources  of 
mismatch  investigated. 
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l.  Array  Signal  Processing 


1.1.  Introduction 

This  chapter  presents  a  comprehensive  discussion  of  array  processing 
techniques  used  to  process  the  data  from  line  arrays  with  equi-spaced  sensors.  The 
advantage  of  using  an  array  of  sensors  is  the  spatial  discrimination  that  it  provides. 
The  directions  from  which  underwater  signals  propagate  can  be  extracted  from  the 
measured  data.  After  an  outline  of  the  frequency-wavenumber  representation  of  the 
acoustic  pressure  field  and  its  associated  assumptions,  the  conventional 
beamformer  is  derived.  Beamforming  is  decomposed  into  a  signal  model  and  a 
processing  structure.  The  signal  model  provides  a  mathematical  description  of  the 
propagation  physics  and  makes  assumptions  on  the  arrival  wavefronts.  The 
function  of  the  processing  structure  is  to  measure  how  close  the  experimental  data 
are  to  a  particular  direction  of  arrival.  The  adaptive  Minimum  Variance  processor  is 
emphasized  and  will  be  used  in  the  thesis.  Its  implementation  is  discussed  and  its 
performance  summarized. 


1.2.  Basic  Assumptions 

Array  processing  methods  decompose  the  pressure  field  into  a  directionality 
representation  by  assuming  homogeneous  wavefields  [Yen;  1977],  A  homogeneous 
wavefield  is  the  superposition  of  uncorrelated  elementary  waves  [Bohme;  1987], 

This  mathematical  formulation  sets  a  clear  analogy  between  time  series  spectral 
estimation  and  array  processing  in  its  frequency-wavenumber  representation  of  the 
random  field  [Bohme;  1987.  Capon;  1969]  : 

p{  t,  r)  =  /  exp[  j  ( -  u  t  +  Fr) )  dZ{u,  F)  (1.1) 

where  p  denotes  the  random  field  at  time  f  and  spatial  position  r,  w  the  circular 
frequency  in  rad/s  and  F  the  wavenumber  vector.  dZ  is  the  spectral  increment  of 
the  random  spectral  measure  Z(w  ,  F)  with  covariance 

COv\dZ{ w  ,  F)dZ(w',  F  )]  =  -  5(w— w'  )  <5(  F—  F)  S(w,F)  du  dF  (1.2) 

t  J  (2jt  r 

where  w  and  w'  are  circular  frequencies,  F  and  F  wavenumber  vectors  ,  S  the 
frequency-wavenumber  spectrum  of  the  random  field  p(t  ,  r )  and  5  the  dirac  delta 
function. 

This  representation  is  a  generalization  of  the  Cramer  spectral  representation 
for  stationary  time  series  [Priestley;  1981].  Its  main  assumption  is  that  spectral 
increments  dZ  at  two  distinct  frequencies  or  at  two  distinct  wavenumber  vectors 
are  uncorrelated.  It  constitutes  the  foundation  of  any  array  processing  method. 
Thus,  this  representation  is  one  where  the  receiving  array  measures  signals  which 
result  from  the  superposition  of  uncorrelated  elementary  waves.  These  elementary 
waves  generally  are  assumed  planar.  Such  a  homogeneous  wavefield  representation 
is  restrictive  since  it  excludes  correlated  arrivals  (e.g.  multipaths). 


1.3.  Conventional  Beamforming 

Conventional  beamforming  on  a  line  array  with  uniformely  spaced  sensors  is  a 
simple  illustration  of  the  already  mentioned  analogy  with  classic  time  series 
spectral  estimation  [Kay;  1988].  The  conventional  beamformer  is  the  simplest 
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processing  structure  since  it  implements  a  linear  regression  on  the  sensor  outputs, 
either  in  the  time  domain  or  in  the  frequency  domain  [Shumway;  1988]. 

The  time  domain  implementation  is  known  as  delay-and-sum  beamforming.  A 
beam  is  formed  in  the  time  domain  by  summing  up  the  lagged  (delayed)  outputs  of 
the  array  sensors  [Johnson;  1982,  De  Fatta;  1988] 

M— 1 

y(t,6)=  £amxm(t-rm)  (1.3) 

m—  0 

where  t  denotes  time,  xm(t)  is  the  output  of  the  mth  sensor,  M  the  number  of 
sensors,  (am)mm0.  m-i  the  shading  coefficients  and  0-m)m_o,  m-i  the  delays  associated 
with  each  sensor  corresponding  to  a  particular  look  direction,  9.  The  time  delays  t„ 
are  determined  by  making  additional  assumptions  on  the  propagating  signals.  If 
one  considers  plane  waves  impinging  on  a  line  array  with  equally  spaced  sensors, 
the  delays  t„  are  given  by 

m  d  sin0  .. 

Tm  - - - -  (1.4) 

where  d  is  the  interelement  spacing,  c  the  sound  velocity  and  9  the  look  direction 
[DeFatta;  1988]. 

The  frequency  domain  regression  can  be  derived  by  taking  the  Fourier 
transform  of  the  output  y(t) 


Y(f,9)=  E  ame 

m—  0 


-j2ir  fr 


Xm(f) 


where  f  is  the  carrier  frequency,  Y(  f,  9)  and  Xm(  f)  the  Fourier  transforms  of  y(r) 
and  xm(t)  respectively.  In  essence,  Equation  (1.5)  is  the  narrow-band  implementation 
of  Equation  (1.3)  at  the  frequency  f,  and  can  also  be  expressed  as  an  inner  product : 

Y(f.k)  =  AH  X  (1.6) 


xr=  [xa(f)  xx(f) 


4-1  (f) 


a  T  j2xfrl  j2xfru 

A  =  \a0e  10  axeJ  1  •••  aM_xeJ  "~l  (1.8) 

the  superscript  H  denotes  the  Hermitian  operation  (complex  conjugate  transpose 
operation),  and  the  superscript  T  denotes  the  transpose  operation.  Assuming  plane 
waves,  that  is  using  Equation  (1.4)  and  letting 

k.-L^SL  (1.9) 


Equation  (1.5)  can  be  put  under  the  form  of  a  spatially  windowed  Discrete  Fourier 
Transform  [Oppenheim;  1975]  : 


Y(f,k)=  £am  e~j2irkm  Xm(f) 

mm  0 


(1.10) 


where  k  is  a  normalized  wavenumber.  The  frequency-wavenumber  power  spectrum 
is  estimated  by  the  periodogram 


lp{f.k)  =  - - 1  Y{f,k)  I2 

(  E  a,)2 
o 


(1.11) 
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where  M  is  the  number  of  sensors.  Ip  (f,  k)  can  be  converted  into  a  power 
spectrum  parametrized  by  the  arrival  angle  9  by  the  nonlinear  mapping  given  in 
Equation  (1.9)  (De  Fatta;  1 9S8I-  The  narrow-band  conventional  beamformer  is  often 
implemented  as  a  FFT  beamformer  (Williams;  1968].  Since  it  is  well  known  that  the 
periodogram  gives  an  inconsistent  estimate  of  the  power  spectral  density  [Priestley; 
1981],  additional  incoherent  averaging  over  successive  time  snapshots  is  typically 
performed.  The  estimate  then  becomes 


P(f.k)  = 


1 


<"£V  K'-° 

i— 0 


i  K-l  2 

1  E  I  Yj{  f.k)  I 


(1-12) 


where  f.k)  is  FFT  beamformer  output  from  the  ith  time  snapshot.  P  has  a  \2 
distribution  with  2  K  degrees  of  freedom  [Kay;  1988].  The  mean  and  the  variance  of 
the  conventional  estimator  are  [Pillai;  1989,  p.  112] 

E(P)  =  P  (1.13a) 

Var(P)  =  jP2  (1.13b) 

where  E  denotes  the  expectation,  Var  the  variance,  and  P  the  power  spectrum 
computed  by  the  conventional  processor  on  the  expectation  of  the  covariance 
matrix. 

The  conventional  FFT  beamformer  directly  operates  on  the  data,  and 
corresponds  to  a  direct  form  of  spectral  estimation.  The  power  spectrum  can  also 
be  obtained  by  indirect  transformation  of  the  cross-spectral  matrix  based  on  the 
data  [Marple;  1987].  This  formulation  of  the  conventional  beamformer  in  terms  of 
the  array  covariance  matrix  is  now  examined. 

When  the  window  applied  to  the  aperture  is  the  rectangular  window,  the  afs 
are  taken  equal  to  unity  and  the  power  spectrum  can  be  expressed  in  a  simplified 
matrix  form 

Ip  (f  .  k  )  =  (EHX)  (EHX)H  (1.14a) 


EhX  XmE 


M2 


Hi 


(1.14b) 


where  E  is  the  plane  wave  steering  vector, 

E7'=[l  e./2*k  ...  gj  2ir  k  (M— 1) 


(1-15) 


and  X  XH  is  an  estimate  at  the  frequency  f  of  the  cross  spectral  density  matrix 
based  on  the  data,  or  array  covariance  matrix.  With  additional  averaging  over  time 
snapshots,  the  conventional  beamformer,  also  called  Bartlett  beamformer,  is  given 
by 

PB(f.k)=~^EHRE  (1.16) 


with 

R  =  -pr  *E*X,  X"  (1.17) 

ft  i  -  o 

where  X,  is  the  vector  of  Fourier  coefficients  corresponding  to  the  time  snapshot  /. 
R  is  an  estimate  of  the  array  covariance  matrix.  When  a  weighting  window  is 
appplied,  one  gets 
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(1.18) 


EH  WH  R  W  E 


where  the  (M,M)  real  diagonal  matrix  W  is 


Wr  = 


1 


i  M— 1 
77  £ 

Mi-  o 


0/ 


a o 

0 

...  o 

0 

a\ 

•  •  • 

• 

• 

...  o 

0 

0 

'  ‘ '  <*m- i 

(1.19) 


1.4.  Matched  Filter  Interpretation  of  Beamforming 

This  indirect  formulation  of  the  conventional  beamformer  in  terms  of  the  array 
covariance  matrix  allows  the  decomposition  of  the  processing  into  two  components: 
a  processor  and  a  set  of  steering  vectors. 

For  a  given  set  of  replica  vectors,  finding  an  arrival  corresponds  to  finding  a 
peak  in  the  processor  output.  Therefore,  the  processing  can  be  interpreted  as 
spatial  matched  filtering.  Using  Equation  (1.14),  the  quantity  to  be  maximized  in 
order  to  find  an  arrival  or  a  signal,  is  given  by 

EH  X  XH  E  =  IEH  XI2  (1.20) 

where  E  is  a  direction  vector  (given  by  Equation  (1.15))  corresponding  to  an  arrival 
angle  0.  E  is  a  steering  vector  or  replica  vector  from  a  set  often  called  array 
manifold  [Schmidt;  1981,  Clarke;  1988]  which  contains  the  true  underlying  arrival 
wavefronts.  Generally,  the  plane  wave  array  manifold  is  selected. 

If  a  single  signal  impinges  on  the  array  at  a  given  but  unknown  angle  of  arrival, 
the  measured  array  signal  vector  X  given  by  Equation  (1.7)  will  be  an  amplitude 
weighted  plane  wavefront  vector.  From  the  Schwartz  inequality  the  output  of  the 
processor  will  be  maximum  if  the  replica  vector  E  and  the  array  signal  vector  X  are 
colinear.  The  processor  provides  a  measure  of  the  mismatch  between  the  observed 
field  and  a  replica  vector  or  steering  vector.  A  match  implies  a  high  output  and  an 
arrival  [Cox;  1973]. 

Most  of  the  emphasis  in  the  array  processing  community  has  been  on  the 
choice  of  a  best  processing  structure  and  has  lead  to  the  development  of  various 
techniques  besides  conventional  beamforming,  such  as  adaptive  methods, 
approaches  based  on  auto-regressive  models,  and  principle  component  techniques 
[Pillai;  1989].  All  these  advanced  processing  schemes  operate  on  the  array 
covariance  matrix  and  have  their  own  advantages  and  disadvantages  related  to 
performance  criteria  such  as  resolution,  processing  complexity,  stability.  In 
particular,  the  adaptive  Minimum  Variance  method  [Capon;  1969]  has  received  much 
attention  in  the  underwater  acoustics  community  and  will  be  studied  later  in  this 
chapter  and  used  throughout  this  thesis.  The  propagation  physics  play  a  minor  role 
in  conventional  beamforming  since  it  is  incorporated  through  a  simple  signal  model 
and  the  selection  of  planar  wavefronts.  This  underlying  signal  model  is  now 
studied. 
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1.5.  Signal  Model  Definition 

If  p  narrowband  signals  at  the  carrier  frequency  f  impinge  on  the  array,  the 
time  domain  array  signal  vector  at  the  input  of  the  beamformer  is  given  by  [Cox; 
1973,  De  Graaf;  1985] 

r(f)  =  B  s(f)  +  n(t)  (1.21) 

where  s (t)  is  the  (p,  1)  signal  vector,  sometimes  called  the  envelope  vector  [Cadzow; 
1988],  n(r)  is  the  (M,l)  vector  of  additive  noise  and  B  the  matrix  of  stacked 
directions  with  columns  equal  to 

B=  [d0  d,  ...  d,_,]  (1.22) 

each  (M.l)  column  vector  d j  corresponds  to  an  arrival  angle  ffj.  The  equivalent 
narrow-band  frequency  domain  representation  is,  at  the  signal  carrier  frequency  f, 

r  =  Bs  +  n  (1.23) 

where  s  is  the  (p,\)  signal  vector,  n  is  the  additive  noise  (M,l)  vector,  and  B  the 
matrix  of  stacked  directions.  The  corresponding  model  covariance  matrix 
R model  -  r  rH  is  thus 

RWe,  =  BSB«  +  Q  (1.24) 

where  S  is  the  ( p,p )  signal  covariance  matrix  given  by  £( s  sH)  and  Q.is  the  noise 
covariance  matrix  given  by  E(n  nH)  with  £  denoting  the  expectation. 

Under  the  assumption  of  a  homogeneous  wavefield,  the  signals  are 
uncorrelated,  and  the  ( p,p )  signal  matrix  S  is  diagonal  of  full  rank  p.  Recalling  that 
S  =  £(s  sH)  with  sr  =  [s0  Si  •  •  •  s„_i],  uncorrelated  arrivals  means  that  the  signal 
envelopes  are  such  that,  for  j, 

E(Sj  sj)  =  0  (1.25) 

where  *  denotes  the  complex  conjugate  operation.  If  s,  and  Sj  are  correlated,  then 

E(Sj  sj)  ^  0  (1.26) 

Restricting  the  discussion  to  two  arrivals,  the  signal  covariance  matrix  is 

E(s0So)  E(s0s\) 

E(sq  Si)  E(si  s] ) 


(1.27) 


with  £(s,5,*)  =  <t },  the  power  of  the  ith  signal,  and  for  /  ^  j,  E{Sj  sj)  =  (TjOjPij,  where 
Pi.j  ~  Pj.i 's  (he  complex  correlation  given  by 


_  £(s<  sj) 

PiJ~  £(ISjl2)£(ISjl2) 


(1.28) 


Then,  the  signal  covariance  matrix  S  can  be  written  [Reddy;  1987] 


For  fully  correlated  arrivals  Ipl  =  1,  hence  det  S  =  <t%  a\  (1  -  \p\2)  =  0.  Fully 
correlated  arrivals  are  called  coherent  arrivals  [Piliai;  1989,  p.  10]. 
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The  standard  model  for  the  noise  covariance  Q.is  a  sensor  noise  uncorrelated 
from  sensor  to  sensor  [Nickel;  1988,  Schmidt;  1986).  Under  this  assumption,  the 
noise  covariance  matrix  is  Q.=  <r;j  I,  where  a2  is  the  noise  power  and  I  the  identity 
matrix.  It  follows  that  the  (M.  M)  covariance  matrix  R model  is  of  full  rank. 

Furthermore,  if  one  assumes  the  elementary  waves  of  the  field  to  be  planar, 
the  covariance  matrix  Rmodei ,  which  is  Hermitian  by  construction  (i.e. 

RH model  =  Rmodei  )•  also  is  Toeplitz  (constant  along  the  diagonals). 

R model  =  ffo2  d0  d0H  +  a I2  d!  djH  +  •  ■  •  +  op_?  dp_!  d,.^  +  a2  I  (1.30) 


Since  d,  =  1  ej2nk'  eJ2*2ki  ■  •  •  eJ2n{M  l)k‘  J,  each  outer  product  in  the  sum  given  by 
Equation  (1.30)  is  Toeplitz  : 


d  dfH  = 

I  i 


1 

J2nki 


e~J2*k‘  e~j2r2k‘ 

1  0-i2l,k< 


J2idM-2)ki 

e  e  ' 


—j2HM—l)kl 


(1.31) 


so  that  R modei  is  both  Toeplitz  and  Hermitian.  The  Toeplitz  character  of  the  array 
covariance  matrix  is  a  result  of  the  duality  between  homogeneous  plane  wavefields 
and  wide  sense  stationary  time  series. 


1.6.  Adaptive  Beamforming 

The  narrow-band  beamforming  problem  can  be  viewed  as  finding  a  set  of 
complex  weights  or,  equivalently,  a  steering  vector  given  by  Equation  (1.15) 

[De  Graaf;  1985]  corresponding  to  the  signal  direction  of  arrival.  The  conventional 
processor  (which  assumes  plane  waves)  is  the  most  widely  used,  the  best  known 
and  the  most  robust,  but  has  the  least  resolution.  Since  only  limited  sampling  in 
space  can  be  achieved  at  a  reasonnable  cost,  an  important  performance  criteria  is 
resolution.  Data  adaptive  techniques,  which  are  described  in  what  follows,  are  used 
to  meet  this  design  requirement.  They  make  full  use  of  the  information  summarized 
in  the  covariance  matrix,  using  the  signal  model  described  in  the  previous 
paragraph. 


1.6.1.  Optimum  versus  Adaptive  Processing 

Optimum  signal  processing  derives  optimal  processors  based  on  completely 
known  covariance  matrices  [Orfanidis,  1988].  It  assumes  that  all  the  entries  in  the 
signal  model  given  by  Equation  (1.24)  are  perfectly  known,  in  particular  the  noise 
covariance  matrix  Q,  Adaptive  arrays  which  find  multiple  applications  in  RADAR  or 
active  SONAR  [Brennan;  1973],  solve  the  problem  of  detecting  a  signal  in  a  "primary 
sample"  based  on  a  signal-free  or  noise-only  "secondary  sample”  [Kelly;  1989].  The 
secondary  samples  allow  the  estimation  of  the  noise  covariance  matrix  Q,  which  is 
used  to  compute  the  optimal  set  of  weights  which  maximizes  the  output  signal-to- 
noise  ratio,  given  a  particular  direction  of  look  E  [Cox;  1973]: 


w  Q--1  E 

ww  —  — — - ; 

eh  or'  e 


(1.32) 


The  processor  using  those  weights  has  been  normalized  to  yield  unity  response 
with  zero  phase  shift  in  the  look  direction  E  [Cox;  1987].  In  adaptive  arrays,  the 
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look  direction  E  is  called  unadapted  beam  or  quiescent  beam.  The  output  power 
produced  by  the  processor  is  given  by 

Pprocnsor  =  WH  R  W  (133) 


In  the  context  of  passive  SONAR,  the  noise  field  includes  all  "unwanted"  signals 
and  the  noise  covariance  matrix  Qis  not  known.  Then,  the  optimum  processor  is 
implemented  by  an  adaptive  processor  which  uses  estimate  of  the  array  covariance 
[Cox;  1987]  (since  the  noise  covariance  matrix  is  not  known,  the  estimate  of  the 
whole  array  covariance  R  is  used  instead  of  the  noise-only  covariance  matrix  Q). 
Then,  the  adaptive  processor  is  obtained  by  using  the  weight 


W  = 


R-1  E 

eh  r-1  E 


(134) 


and  is  the  well  known  Minimum  Variance  Distorsionless  Response  (MVDR) 
beamformer  [Zoltowski;  1988].  This  adaptive  processor  is  called  the  Minimum 
Variance  (MV)  processor  in  this  thesis.  It  was  originally  called  Maximum  Likelihood 
Method  (MLM)  [Capon;  1969]  because  it  was  shown  to  be  the  maximum  likelihood 
estimator  in  the  case  of  Gaussian  noise  (i.e.  it  minimizes  the  log-likelihood  under 
the  assumption  of  Gaussian  statistics).  The  MV  processor  is  also  called  the 
Minimum  Energy  processor  in  [De  Graaf;  1985].  This  adaptive  beamformer 
determines  a  data  dependent  set  of  weights  as  the  solution  of  a  constrained 
optimization  problem  [Griffiths;  1977,  Monzingo;  1985,  Van  Veen;  1988].  The  MV 
processor,  classically,  is  derived  by  minimizing  the  output  power  of  the  adaptive 
spatial  filter  subject  to  the  constraint  of  passing  the  on-look  direction  undistorted 
[Lacoss;  1971,  Cox;  1973,  Kanasewich;  1975,  Marple;  1987]  : 


MIN 


\k/H  n  ui 
«  K  W  J£hw  _  j 


(1.35) 


where  W  is  the  unknown  weight  vector,  E  the  look  of  direction  vector  or  replica 
vector  from  the  array  manifold  and  R  the  cross-spectral  density  matrix,  estimated 
from  the  data  and  given  by  Equation  (1.17).  As  before,  one  assumes  that  the  data 
have  been  Fourier  analysed  and  a  single  frequency  f  is  considered.  The  weight 
vector  W  can  be  derived  by  using  the  Lagrange’s  multiplier  method  [Lacoss;  1971, 
Kanasewitch;  1975,  Marple;  1987]  and  is  given  by  Equation  (1.34).  The 
corresponding  power  spectrum  is 

Pm(  f.  k)  =  -eh^-ie  (1-36) 


As  stated  earlier,  this  processor  also  maximizes  the  output  signal-to-noise  ratio 
and  the  array  gain  (the  output  signal-to-noise  ratio  in  the  beam  divided  by  the 
signal-to-noise  ratio  at  a  sensor)  [Gammelsaeter;  1980].  If 

R  =  crs2d  dH+Q  (1.37) 


where  as2  is  the  signal  power  and  d  the  signal  direction  vector  (normalized  to  unit 
norm),  the  array  gain  is  given  by 


(Tj2  WH  d  dH  W 

/ 

°s2 

wH  aw 

/ 

(138) 


where  qx  ]  is  the  (1,1)  term  of  the  noise  covariance  matrix  Q,  Then, 

ql  !  (wH  d  dH  w] 

°=  WhQ.W 


(1.39) 
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which  can  be  transformed  into 


G 


(»T  i 


WH  d  dH  W 

WH  R  W  -  a2  WHddHW 


(1.40) 


It  is  evident  from  Equations  1.39  and  1.40  that  maximizing  the  array  gain  G  is 
equivalent  to  minimizing  WH  QW  or  minimizing  WH  R  W  if  one  constrains  WH  d  to 
unity. 

Assuming  the  signal  model  of  Equation  (1.37)  and  spatially  white  noise  i.e. 

Q.=  a2 1,  the  MV  power  spectrum  can  be  expressed  as 


Pmv  = 


M 


*„2(1  +  % 


1  +  \EH  dl: 

a2  M 


(1.41) 


using  the  matrix  inversion  lemma  [Kay;  1988,  p.  24].  If  there  is  no  mismatch.  E 

corresponds  to  the  signal  direction  and  IEH  dl  =  1  so  that  PMV  =  -^(<r|  +  <r2).  It  is 

identical  to  the  output  of  the  conventional  processor  under  the  same  assumptions 

since  PB  =  -tt-Eh  R  E  =  IEH  dl2  -f  <r%].  Similarly,  in  a  noise  direction  EH  d  is 

M  M  M 

small  so  PMV  and  PB  are  close  to  -7-  a2. 

M 

The  MV  processor  stands  as  a  "high  resolution"  method  compared  to  the 
conventional  FFT  beamformer.  This  statement  is  proven  by  showing  that  PMV  <  PB 
[Pillai;  1989].  The  covariance  matrix  R  is  Hermitian  (RH  =  R),  and  can  be  decomposed 
according  to  the  spectral  theorem  [Golub;  1985]  as 

R  =  aHDQ.  (1.42) 

where  Qis  unitary  (i.e.  Q.H  Q.  =  I,  where  I  is  the  identity  matrix),  and  D  is  a  diagonal 
real  matrix.  The  covariance  matrix  can  also  be  put  under  the  form 

R  =  (D*Q)hD*Q.=  (R*)hR*  (1.43) 

Assume  the  steering  vector  e  normalized  such  that  eH  e  =  1,  the  unit  norm  relation 
can  be  expressed  as 


eHR*R-*e=l  (1.44) 

which  is  a  vector  inner-product  between  two  vectors  u  =  (R*)H  e  and  v  =  R~Kj  e. 
Using  the  Schwarz  inequality  uHv  <  uHu  vHv  and 

uHu  * eHRe 


vHv  = 

one  obtains  the  inequality  [Pillai;  1989,  p.  20] 

1 


e«R-‘-  - 


<  eHRe 


(1.45a) 

(1.45b) 

(1.46) 


so  that,  while  the  MV  processor  reports  the  same  peak  levels  as  the  conventional 
Bartlett  processor,  it  has  a  better  resolving  power  in  general.  That  is  why  the  MV 
estimator  is  called  a  "high  resolution"  method  and  is  used  in  various  areas  such  as 
the  seismic  field  [Capon;  1969,  Baggeroer;  1974,  Baggeroer;  1982,  Hsu;  1986],  the 
oceanographic  field  [Oltman-shay;  1984]  as  well  as  the  underwater  acoustics  arena 
[Owsley;  1985]. 
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1.6.2.  Behaviour  of  the  MV  Beamformer 

The  behaviour  of  the  MV  beamformer  has  been  studied  rather  extensively  [Cox; 
1973,  De  Graaf;  1985]  and  only  general  performance  trends  are  presented  here.  For 
the  adaptive  array  pointing  in  a  given  direction,  the  signal  arriving  in  another 
direction  is  an  interference  to  be  cancelled.  The  stronger  the  interference,  the  more 
suppression  the  designed  filter  is  able  to  put  on  it.  This  suppression,  based  on  the 
covariance  structure  across  the  array,  depends  on  the  signal-to-noise  ratio.  In  the 
context  of  matched  filtering,  the  MV  processor  is  such  that  a  slight  mismatch  in  the 
replica  vector  results  in  serious  reduction  of  its  output,  in  .he  case  of  a  high 
signal-to-noise  ratio  environment  [Cox;  1973).  This  high  sensitivity  to  mismatch  can 
cause  problems,  for  instance,  when  an  insufficient  number  of  look  directions  is 
used  or  if  the  array  manifold  is  chosen  with  improper  assumptions  about  the  true 
underlying  wavefronts.  Sensitivity  to  mismatch  also  means  that  any  kind  of  errors 
(such  as  gain  or  phase  errors  in  the  analog  processing  chain  of  each  sensor  or  such 
as  sensor  position)  will  degrade  the  performance  of  the  adaptive  processor  to  a 
large  extend.  Then,  the  MV  beamformer  breaks  down  and  its  resolving  capability 
can  be  less  than  that  of  the  conventional  processor  [Seglison;  1970,  McDonough; 
1972],  As  a  result  of  this  sensitivity  to  mismatch,  resolution  is  highly  dependent  on 
signal-to-noise  ratio,  and  is  high  for  high  signal-to-noise  ratios.  For  instance,  the  MV 
resolution  for  a  line  array  receiving  two  signals  with  30  dB  output  signal-to-noise 
ratio  is  slightly  over  three  times  the  Rayleigh  resolution  of  conventional  processor 
[Cox;  1973]  which  is  given  in  radians  by  [Burdic;  1984] 

0  =  j  (1.47) 

where  X  is  the  acoustic  wavelength  and  L  is  the  length  of  the  array.  At  low  signal- 
to-noise  ratio,  the  sensitivity  to  mismatch  of  the  MV  processor  decreases  and  the 
resolution  of  the  adaptive  processor  becomes  comparable  to  the  conventional 
processor  resolution  [Cox;  1973,  Kanasewich;  1975,  Zoltowski;  1988]. 

1.6.3.  Implementation  Issues 

The  MV  beamformer  is  a  more  complicated  processing  structure  than  the 
conventional  processor  and,  to  obtain  the  power  spectrum,  requires  the  estimation 
of  the  narrowband  covariance  matrix  based  on  the  data  and  its  inversion 

1. 6.3.1.  Estimating  the  Covariance  Matrix 

The  problem  of  estimating  the  covariance  matrix  is  avoided  by  the 
conventional  FFT  beamformer.  Since  the  MV  processor  requires  the  inversion  of  the 
covariance  matrix,  a  requirement  is  to  work  with  a  nonsingular  covariance  matrix. 
Thus,  in  Equation  (1.17),  one  must  average  at  least  as  many  dyad  products  as 
sensors,  otherwise  the  covariance  matrix  theoretically  is  singular  and  non- 
invertible.  The  way  one  estimates  the  covariance  matrix  has  some  impact  on  the 
final  power  spectrum. 

This  issue  has  been  addressed  by  the  adaptive  array  community,  and  the 
literature  provides  guidelines  on  how  the  covariance  matrix  should  be  estimated. 
These  guidelines,  derived  for  the  noise  covariance  matrix  CL  immediately  apply  for 
the  full  array  covariance  matrix  R,  since  in  the  context  of  passive  SONAR,  noise 
cannot  be  differentiated  from  signals.  As  already  stated  by  Equation  (1.17),  the 
noise  covariance  matrix  is  estimated  by  averaging  dyad  vectors 

Q=  X.H  0-48) 

K  i  -  o 
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where  X;  is  a  complex  sample  vector  of  the  M  sensor  outputs  corresponding  to  the 
frequency  f.  K  is  the  number  of  snapshots  such  that  K  >  M.  It  is  assumed  that  the 
sensor  outputs  X/  are  independent  multidimensional  complex  Gaussian  random 
vectors  with  probability  density  function  given  by 

fx(*i)  =  7ditaexp  ( ~  Xt"  a_1  X,1  (L49) 

where  the  noise  covariance  matrix  Q.is  Hermitian  positive  definite  and  detQ.is  the 
determinant  of  Q,  It  is  well  known  [Goodman;  1963]  that  the  estimator  given  in 
Equation  (1.17)  or  (1.48)  is  the  maximum  likelihood  estimator  based  on  the  data. 

That  is,  it  maximizes  the  log  likelihood  function,  In  L,  where  L  is  the  joint 
probability  of  the  K  independent  random  vectors  (X,),  _0,  k  -  i-  The  distribution  of 
the  covariance  matrix  estimate  is  a  complex  Wishart  distribution  [Goodman;  1963]. 
Based  on  this  result,  it  has  been  shown  [Reed;  1974]  that  the  adaptive  signal-to- 
interference-plus-noise  ratio  achieved  with  the  weights  W  =  k  Q"1  E  (where  k  is  a 
complex  number,  in  our  discussion  k  =  1/(EH  Q-1  E) )  approaches  the  optimum  one 
(where  Q.is  assumed  known)  by  a  normalized  signal-to-noise  ratio  p  with  expected 
value  given  by 

E  (p)  =  (1.50) 

Equation  (1.50)  leads  to  the  useful  rule  of  thumb  that  one  must  average  a  number  of 
dyads  equal  to  twice  the  number  of  sensors  to  compute  the  covariance  matrix  in 
order  to  maintain  an  average  loss  ratio  of  3  dB.  The  analysis  is  developed  further  in 
[Boroson;  1980]  where  it  is  pointed  out  that  a  larger  number  of  averages  is  required 
to  insure  a  certain  probability  of  having  a  reduction  in  signal-to-noise  ratio  below  3 
dB  (the  examples  given  in  [Boroson;  1980]  are  K  >  3  M  to  ensure  that  the 
probability  P(p  >  3  dB)  <  0.0196,  and  K  >  4  M  to  ensure  that  the  probability 
P(p  >3  dB)<  0.0032). 

It  noted  in  [Gabriel;  1980,  Hudson;  1980]  that  averaging  over  only  a  few 
snapshots,  still,  can  lead  to  an  adaptive  filter  with  good  interference  rejection 
capabilities,  but  with  a  distorted  beampattern.  The  asymptotic  sidelobe  level  of  the 

adaptive  filter  is  given  by  *  [Kelly;  1989].  If  the  noise  plus  interference  field  is 

x  +  1 

composed  of  only  a  few  interferences  with  power  levels  well  above  the  noise 
background,  the  covariance  matrix  R  will  have  a  large  spread  in  eigenvalues  an  a 
therefore  will  be  poorly  conditioned  and,  in  fact,  intrisincally  singular  [Hudson; 

1981,  p.  136].  In  this  case,  the  amount  of  averaging  required  to  estimate  the 
covariance  matrix  can  be  relaxed  to  twice  the  number  of  strong  interferences 
[Hudson;  1981,  p.  123].  Even  with  few  averages,  virtually  perfect  nulling  can  be 
achieved:  the  sample  covariance  matrix  is  rendered  invertible  by  augmenting  its 
diagonal  by  a  supplement  e  [Capon;  1969,  Hudson;  1981,  Carlson;  1988]. 

The  number  of  snapshots  entering  in  the  estimation  of  the  array  covarian.e 
matrix  controls  how  well  the  covariance  matrix  is  estimated,  thus  the  bias  and 
variance  of  the  MV  estimator.  Based  on  the  distribution  of  the  covariance  matrix 
estimate,  the  MV  power  spectrum  has  been  shown  in  [Capon;  1970]  to  have  a  \2 
distribution  with  a  number  of  degrees  of  freedom  equal  to  2(/C-M+l)  where  K  is 
the  number  of  snapshots  and  M  the  number  of  sensors.  The  MV  power  spectrum 
estimator  has  a  mean  and  a  variance  given  by  [Pillai;  1989,  p.  112] 

E{PMv)  =  K~~k  +  1  p  d-Sla) 

Var(PMV)  =  K  -  ^  +  1  P2  (1.51b) 
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where  P  the  power  spectrum  computed  by  the  MV  processor  on  the  expectation  of 
the  covariance  matrix.  Therefore,  the  MV  estimator  is  biased  by 

bias(dB)  =  10  log  K  ~  +  -1  (1.52) 

Furthermore,  one  observes  that  the  adaptive  processor  requires  a  large  number  of 
snapshots  compared  to  the  conventional  processor  in  order  to  achieve  statistical 
stability,  especially  for  large  number  of  sensors. 

1.6.3. 2.  Inverting  the  Covariance  Matrix 

The  possibility  of  handling  an  ill-conditioned  covariance  matrix  was  mentioned 
in  the  previous  Section.  This  problem  is  now  addressed.  Two  treatments  of  singular 
or  nearly  singular  matrices  are  proposed  in  the  literature.  One  uses  the  pseudo¬ 
inverse  instead  of  the  true  inverse  (which  in  the  worst  case  does  not  exist)  and  the 
other  stabilizes  the  covariance  matrix. 


Using  the  Pseudo-Inverse 

The  pseudo-inverse  (also  called  Moore-Penrose  inverse  [Penrose;  1985])  of  a  M 
sensor  array  covariance  matrix  R,  with  k  <  M  non  zero  eigenvalues,  is  defined  as 

R+=  EX-le„e»  (1.53) 

n  —  1 

where  X„  and  e„  are  the  eigenvalues  and  eigenvectors  of  R.  If  R  is  of  full  rank  (i.e. 
k  =  M)  then,  the  pseudo-inverse  corresponds  to  the  true  inverse  R+  =  R~‘.  The  use 
of  the  pseudo-inverse  is  considered  in  (Hudson;  1981,  pp  142-143]  for  the  case 
where  only  k  dyad  products  are  used  to  estimate  the  covariance  matrix  in  an 
adaptive  array  context.  Since  W  is  proportional  to  R-1  E,  the  weight  vector  belongs 
to  the  span  of  eigenvectors  of  R  with  smallest  eigenvalues  which  tends  to  be 
orthogonal  to  the  span  of  eigenvectors  of  R  with  largest  eigenvalues  (those 
correspond  to  the  interferences).  If  X[  >  X2  >  •  •  •  XM,  then  the  dominant  term  in  the 
pseudo-inverse  will  be  X^'e,*  eM  H,  XmLiCM_!  eM_,  H,  •  •  •  ,  and  so  on. 

It  is  shown  in  [Hung;  1983]  through  theoretical  derivations  and  simulations 
that  using  the  pseudo-inverse  can  lead  to  misleading  results.  The  theoretical  reason 
is  that  for  a  matrix  whose  smallest  eigenvalues  are  very  close  to  zero,  the  pseudo¬ 
inverse  (which  in  practice  does  not  keep  those  near  zero  eigenvalues)  actually  does 
not  keep  the  dominant  terms  which  appear  in  the  inverse  of  the  covariance  matrix. 
Thus  the  weight  vector  W  is  chosen  in  a  vector  space  which  does  not  include  the 
dominant  null-eigenvalue  vector  space. 


Using  White  Noise  Stabilization 

In  contrast,  the  white  noise  stabilization  procedure  retains  those  dominant 
terms  which  enter  into  the  true  inverse  R_I, 


k  M 

R  +  €  I  =  E  (X,-  +  c)  e,  e,H  +  E  £  e,  e/* 

i  -  1  i -k  +  1 

(R  +  «  I)"1  =  E  (X,  +  ()''  e,  e"  +  S  e"1  e,  e,H 

i  -  1  i  -  k  +  1 


(1.54) 

(1.55) 


so  that  even  if  a  slightly  incorrect  weighting  of  the  various  directions  is  performed, 
the  optimum  weight  vector  W  is  drawn  from  the  subspace  of  the  null  eigenvalues  of 
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R  [Hudson;  1981,  p.  141]. 

The  idea  of  stabilizing  the  covariance  matrix  was  introduced  early  on  by  Capon 
[Capon;  1969].  The  condition  number  of  the  matrix,  which  is  defined  as 


V - / 

Amax 

where  \min  and  Xmax  are  respectively  the  minimum  and  maximum  eigenvalue  of  the 
array  covariance  matrix,  provides  a  convenient  criterion  for  the  matriv  'vertibility 
(based  on  the  eigenvalue  spread).  The  matrix  can  be  easily  inverted  i 

c  >  c0  (1.57) 


where  c0  is  typically  on  the  order  of  10-6. 

Several  stabilization  schemes  have  been  proposed.  One  employs  the 
knowledge  of  the  extrema  eigenvalues  of  R  [Lunde;  1976]  and  uses  as  stabilization, 
if  c  <  c0: 

£  =  ^max  Cq  Xmin  (1-58) 


if  c  >  c0,  then  no  stabilization  is  needed  and  t  =  0.  Another  scheme,  which  requires 
less  computation  and  is  more  commonly  used,  is  to  add  to  the  main  diagonal  of  R 


i  =  7 


tr(R) 

M 


(1.59) 


where  7  is  a  free  parameter  called  fraction  of  noise  with  typical  values  between  10~3 
and  10-1,  and  tr  is  the  trace  operator.  This  procedure  of  stabilizing  the  covariance 
matrix  with  spatially  white  noise  is  also  called  the  "pseudonoise  addition  technique" 
[Gabriel;  1980,  Takao;  1986],  or  "diagonal  loading"  [Carlson;  1988].  This  stabilization 
scheme  will  be  used  in  this  thesis.  In  addition  to  allowing  or  facilitating  the  matrix 
inversion,  diagonal  loading  can  be  shown  to  increase  the  robustness  of  the  adaptive 
processor  although  there  is  a  trade-off  since  stabilization  is  equivalent  to  an 
increase  of  the  noise  power.  Using  a  stabilization  factor  7  of  0.001  corresponds  to 
intoducing  in  the  system  an  uncorrelated  sensor  noise  30  dB  below  the  average 
sensor  power  level.  Thus,  the  MV  angular  spectrum  is  biased  by  the  amount 

Biasm  -  10  log  ".60, 

Furthermore,  the  general  statistics  of  the  stabilized  covariance  matrix  is  not  known, 
nor  is  the  statistics  of  the  corresponding  MV  estimator  [Baggeroer;  1988].  The 
trade-off  associated  with  the  covariance  matrix  stabilization  has  been  studied  in 
[Carlson;  1988].  By  adding  noise  in  the  system,  one  effectively  reduces  the  adaptive 
filter  nulling  capabilities  against  the  weak  interferences  which  have  small 
eigenvalues.  While  excessive  loading  can  result  in  unacceptable  performance 
degradation,  a  low  level  of  stabilization  is  known  to  greatly  improve  the 
performance  of  adaptive  arrays.  The  various  interpretations  of  the  stabilization 
procedure,  given  in  the  literature  to  explain  the  performance  improvement,  are  now 
presented. 

If  one  uses  R  =  R  +  *  I  (where  «  is  a  free  parameter),  instead  of  R  in  the 
derivation  of  the  MV  processor,  one  does  actually  maximize  with  respect  to  W  the 
following  Lagrangian  functional 

A(W)  =  WH  R  W  +  X(WH  E  -  1)  (1.61) 

which,  in  terms  of  R  can  be  expressed  as 
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A(W)  =  WH  R  W  +  X(WH  E  -  1)  +  £  WH  W 


(1.62) 


where  X  is  a  Lagrange  multiplier  [Gabriel;  1980,  Takao;  1986].  e  which  depends  on 
the  average  power  across  the  array  (c.f.  Equation  (1.59))  is  independent  of  the  look 
direction  and  is  not  a  true  Lagrange  multiplier.  Then,  the  stabilization  procedure  is 
interpreted  as  adding  a  "constraint"  on  the  weight  vector  norm.  Adding  a  constraint 
on  the  adaptive  weight  vector  norm  has  been  proposed  to  increase  the  robustness 
of  the  MV  processor  [Maksym;  1979].  Large  weights  are  characteristic  of 
superdirectivity  conditions  where  an  interference  is  close  to  the  look  direction. 

This  robust  version  of  the  MV  is  also  called  the  Bounded  Weight  Vector  Norm 
Method  (BNX)  [Verlardo;  1989].  Note  that  stabilizing  the  covariance  matrix  and 
performing  the  Minimum  Variance  calculation  is  different  from  using  the  BNX 
processor  for  which  in  Equation  (1.62),  both  X  and  e  would  be  Lagrange  multipliers. 

The  additional  "constraint"  on  the  adaptive  weights,  which  corresponds  to  the 
covariance  matrix  stabilization,  can  be  also  transposed  in  the  dual  optimization 
problem  that  the  MV  processor  solves.  This  dual  problem  is  the  maximization  of  the 
array  gain  [Gammelsaeter;  1980]  and  the  stabilization  "constraint"  is  set  on  the 
white  noise  gain  (Cox;  1987].  The  white  noise  gain  is  obtained  by  setting  Q.=  <r2  I  in 
Equation  (1.37) 


Gw  = 


IWH  dl2 

WH  w 


(1.63) 


where  d  is  the  signal  direction  defined  in  Equation  (1.37).  Although  there  is  not  a 
simple  direct  relationship  between  the  amount  of  stabilization  £  and  the  constraint 
on  the  white  noise  gain  Gw  [Cox;  1987],  this  provides  an  interesting  interpretation 
of  the  covariance  matrix  stabilization  in  terms  of  processor  robustness.  Mismatch 
can  be  expressed  as  a  random  perturbation  in  the  covariance  matrix.  Denoting  the 
perturbation  matrix  by  r?  P,  the  signal  model  of  Equation  (1.37)  becomes 

R  =  <r,2ddH  +  T?P  +  Q.  (1.64) 


and  the  robustness  of  the  processor  can  be  measured  by  the  fractional  sensitivity 
of  the  array  gain  to  the  mismatch  or  the  perturbation  matrix  r?  P.  Defining  the 
fractional  sensitivity  of  the  array  gain  as  [Cox;  1987]: 


r  (dC/dn) 
G 


(1.65) 


one  can  show  that 


WHPW 

IWH  dl2 


(1.66) 


which  is  precisely  the  reciprocal  of  the  array  gain  against  noise  with  covariance 
matrix  P.  When  P  =  I,  S  is  the  reciprocal  of  the  white  noise  gain  given  in  Equation 
(1.63).  This  shows  that  diagonal  loading  actually  increases  the  robustness  of  the  MV 
processor  to  errors  which  are  uncorrelated  from  sensor  to  sensor. 


1.7.  Conclusions 

In  this  chapter,  the  general  concepts  of  array  processing  have  been  outlined. 
Beamforming  has  been  decomposed  into  two  parts:  a  processing  structure  and  a 
signal  model.  The  signal  model  provides  the  general  framework  that  characterizes 
the  physical  environment.  The  processing  structure  is  used  to  decide  whether  or 
not  a  particular  arrival  is  received  by  the  line  array.  The  adaptive  Minimum  Variance 
processor  which  can  provide  more  resolution  than  the  conventional  beamformer  is 
emphasized  and  will  be  used  in  the  rest  of  the  dissertation. 
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2.  Beamforming  in  a  Realistic  Environment 


2.1.  Introduction 

In  Chapter  1.  a  simple  signal  model  was  introduced  where  the  pressure  field  is 
the  superposition  of  uncorrelated  plane  wave  arrivals.  This  chapter  deals  with  a 
more  realistic  and  more  complex  description  of  the  underwater  medium  sampled  by 
a  vertical  line  array.  The  effects  of  the  water  column  boundaries  and  the  variations 
of  the  index  of  refraction  with  depth  are  now  studied.  The  surface  and  bottom 
boundaries  produce  multipath  propagation  or  correlated  arrivals  across  the 
receiving  array.  The  resulting  covariance  matrix  is  not  Toeplitz  as  in  the  original 
signal  model  assuming  planar  wavefronts.  The  MV  processor  is  well  known  to  fail 
in  this  scenario  and  requires  some  pre-processing,  which,  in  order  to  remove  spatial 
correlation,  consists  in  averaging  over  space.  Spatial  smoothing  is  a  necessary 
processing  step  before  using  the  adaptive  processor,  although  it  reduces  the 
aperture  length  (thus  the  resolution)  and  increases  the  computational  burden. 

Then,  the  effect  of  a  variable  sound  speed  with  depth  is  examined.  The  sound 
speed  profile  produces  a  waveguide  propagation,  where  the  arrivals  impinging  on  a 
vertical  array  are  not  planar  but  curved.  The  effect  of  wavefront  curvature  is 
studied  using  a  ray  geometric  approach.  The  plane  wavefront  beamformers  under 
curved  data  can  suffer  large  performance  degradation  for  arrivals  with  low  arrival 
angle  with  respect  to  the  horizontal  and  curved  wavefront  beamforming  is 
proposed.  The  distortion  caused  by  the  smoothing  techniques  in  the  context  of  the 
MV  processor  is  addressed.  Mismatch  can  be  avoided  by  using  curved  steering 
vectors  modified  to  account  for  the  effect  of  spatial  smoothing  on  R. 


2.2.  The  Problem  of  Correlated  Arrivals 

2.2.1.  Limitations  of  the  Adaptive  Methods 

Although  array  processing  and  spectral  analysis  on  time  series  are  related,  they 
are  not  identical.  The  surface  and  bottom  boundaries  produce  multipath 
propagation  [Yen;  1977],  and  the  sound  projected  by  a  source  arrives  at  a  vertical 
line  array  through  multipaths  [Urick;  1983].  The  severe  performance  limitations  of 
the  MV  beamformer  under  correlated  arrivals  are  well  known  [White;  1979,  Cantoni; 
1980].  The  adaptive  filter  uses  the  non-look  direction  signals  to  minimize  its  output 
power  and  since  the  signals  are  correlated  with  the  look  direction  signal,  the  latter 
is  suppressed.  This  process  is  the  so-called  signal  cancellation  phenomena  in 
adaptive  antennas  [Widrow;  1982]. 

2.2.2.  Structure  of  the  Covariance  Matrix  Under  Correlated  Arrivals 

In  Section  1.5  ,  it  was  shown  that  coherence  between  two  arrivals  results  in  a 
rank  deficient  signal  covariance  matrix.  This  explains  why,  on  an  algorithmic  level, 
the  adaptive  method  fails.  Since  the  signal  model  presented  in  Chapter  1  breaks 
down,  the  structure  of  the  covariance  matrix  for  inhomogeneous  wavefields  is  now 
investigated. 

The  array  covariance  matrix  of  a  homogeneous  field  using  the  plane  wave 
array  manifold  is  both  Toeplitz  and  Hermitian.  In  a  realistic  environment,  the  field 
is  not  homogeneous  i.e.  not  stationnary  in  space  because  of  multipath  arrivals  and 
the  covariance  matrix  of  an  equi-spaced  sensor  line  array  is  not  Toeplitz  [Gabriel; 
1980]  (it  is  Hermitian  by  construction).  The  model  covariance  matrix  for  two 
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correlated  arrivals  in  spatially  white  sensor  noise,  given  by  Equation  (1.30)  (where 
the  signal  covariance  matrix  is  given  by  Equation  (1.29))  has  terms  along  the  upper 
diagonal  given  by 


R modAU  +  1)  =  4  e  J2*k°  +  <r?  e  J2*ki  + 


b  e. 


j2xi{k0  —  k  j)  j  e~j2irk0 


+  [p '  ej2”i(kl  -ko>]  e~j2irk' 


(2.1) 


for  i  =  1,  M- 1.  The  matrix  is  not  Toeplitz  because  of  the  last  term  due  to 
correlation  (p  ^  0).  Correlation  introduces  a  modulation  term  in  the  diagonals  of  the 
array  covariance  matrix,  as  can  be  observed  on  Figure  2.1  where  plotted  are  the 
magnitude  and  phase  of  the  upper  diagonals  of  the  array  covariance  matrix  for  a 
field  corresponding  to  two  uncorrelated  and  coherent  unit-power  plane  wave 
arrivals.  The  order  of  the  diagonals  can  be  identified  by  their  length.  The  line  array 

(10y  wavelength  long)  has  32  sensors  with  25  m  interelement  spacing.  The  arrivals 

have  -10°  and  -20°  incidence  angles.  The  operating  frequency  is  20  Hz  and  the 
sound  speed  is  1498  m/s. 

Observing  this  non-Toeplitz  character  is  important  since  it  leads,  in  a  natural 
way,  to  the  processing  methods  designed  to  limit  the  negative  effects  of 
correlation:  spatial  smoothing. 


2.2.3.  Pre-Processing  Methods 

The  methods  which  limit  the  effects  of  spatial  correlation  operate  on  the 
covariance  matrix  to  make  it  "more  Toeplitz"  by  performing  an  averaging  operation 
[Linebarger;  1988].  Based  on  the  concept  of  redudancy  [Linebarger;  1988],  the 
straight  arithmetic  and  geometric  averaging  of  the  covariance  matrix  diagonals, 
have  been  proposed  in  the  literature  [Naidu;  1988,  Hsu;  1986].  However,  forcing  a 
Toeplitz  structure  does  not  systematically  result  in  a  non-negative  definite  cross- 
spectral  matrix  and  is  considered  not  viable  [Linebarger;  1988].  The  major 
successful  averaging  techniques  proposed  in  the  literature  are 

(1)  Frequency  averaging  which  is  averaging  covariance  matrices  at  different 
frequencies  [Wang;  1985,  Zhu;  1988].  This  method  is  relevant  only  to 
broadband  situations  and  will  not  be  discussed  here. 

(2)  Spatial  averaging  or  spatial  smoothing,  developed  for  narrowband  problems 
(and  for  a  plane  wavefront  array  manifold)  [Gabriel;  1980,  Shan;  1985,  Su;  1986, 
Reddy;  1987,  Takao;  1987,  Linebarger;  1988].  This  method  is  particularly 
relevant  to  the  problems  discussed  here,  and  will  be  used  in  this  work. 


2.2.3. 1.  Spatial  Smoothing 

The  concept  of  spatial  smoothing  comes  from  the  idea  that  the  relative  phase 
between  two  correlated  arrivals  changes  in  space  [Gabriel;  1980,  Widrow;  1982].  In 
practice,  moving  the  array  in  space  may  not  be  practical.  Spatial  smoothing  is 
achieved  by  averaging  covariance  matrices  estimated  on  subarrays  extracted  from 
the  full  aperture  and  corresponding  to  slightly  different  spatial  position.  It  induces 
a  reduction  of  the  terms  due  to  correlation  [Linebarger;  1988]  because  the  phase 
variations  are  eventually  averaged  to  zero.  The  technique  uses  the  special 
geometry  of  line  arrays  with  equi-spaced  sensors  [Linebarger;  1988]  and  computes 
the  average  of  the  covariance  matrix  estimated  on  subarrays,  where  each  subarray 
shares  all  but  one  of  its  sensors  with  an  adjacent  subarray  [Shan;  1985].  If  the  full 
aperture  array  has  M  sensors,  and  the  subarrays  have  s  sensors,  this  is  equivalent 
to  averaging  M  -  s  +  1  blocks  of  5  by  s  element  extracted  from  the  full  array  M  by 
M  covariance  matrix  along  the  main  diagonal  [Takao;  1987].  It  can  also  be  viewed  as 
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a  low-pass  filtering  operation  [Linebarger;  1988].  Furthermore  averaging  decreases 
the  variance  of  the  estimate  of  the  cross-spectral  matrix. 

Spatial  smoothing  is  implemented  as  follows  [Shan;  1985].  If  X  if  the  full 
aperture  (M,  1)  array  vector,  and  X,  the  ith  (s.l)  subaperture  signal  vector 

X,r  =  [*(/•)  XM(f)  •••  Xi+S(f)]  (2.2) 

the  spatially  smoothed  covariance  matrix  is 

Ms 

R smooth*  =  ~  T.SX,  X,H  (2.3) 

M  —  S  +  1  /  »  0 


To  illustrate  the  decorrelation  process  achieved  by  spatial  smoothing,  the  diagonals 
of  the  array  covariance  of  a  field  composed  of  two  coherent  arrivals  are  plotted  in 
Figure  2.2  for  two  smoothing  scenarios.  As  in  Figure  2.1,  the  two  unit-power  plane 
wave  arrivals  impinge  on  the  array  with  incidence  angles  of  -  10°  and  -  20°.  The 

array  is  10  y  wavelength  long  and  has  32  sensors.  Two  different  subarray  lengths 

are  selected,  24  and  16,  so  that  the  number  of  averages  in  the  smoothing  increases 
from  9  to  17.  As  the  number  of  average  increases  or  subarray  length  decreases,  the 
covariance  matrix  becomes  more  and  more  Toeplitz. 

It  is  shown  in  [Shan;  1985]  that  if  p  coherent  signals  impinge  on  the  array,  at 
least  p  averages  over  subarrays  of  p  elements  must  be  performed  in  order  to 
recover  the  rank  of  the  signal  covariance  matrix,  that  is  to  have  det  S  *  0  (i.e.  S  is  a 
p  by  p  matrix  of  rank  p).  This  result  is  of  great  importance  because  it  indicates  the 
number  of  averages  necessary  to  determine  the  direction  of  arrivals  with  high 
resolution  methods  (i.e.  eigenvector  methods).  It  leads  to  the  rule  of  thumb  that  a  M 

sensor  array  is  able  to  detect  at  most  —  correlated  arrivals  with  spatial  smoothing 

[Shan;  1985]. 

This  results  is  not  as  valuable  for  the  MV  processor  since  a  full  rank  matrix 
does  not  ensure  total  decorrelation  of  the  arrivals  (i.e.  a  diagonal  signal  covariance 
matrix).  Partial  correlation  among  the  arrivals  leads  to  partial  signal  cancellation  for 
the  adaptive  beamformer  and  therefore  performance  degradation. 

The  spatial  smoothing  decorrelation  rate  for  two  correlated  arrivals  is  studied 
in  [Reddy;  1987].  The  model  covariance  matrix  corresponding  to  the  kth  subarray 
can  be  expressed  as 

R*  =  B  S  (<J>i-I)H  BH  +  <r2  I  (2.4) 


where  Rk  is  the  subarray  covariance  matrix  and  corresponds  to  Xk  Xtf,  B  is  the 
stacked  plane  wave  direction  matrix  (a  s  by  p  matrix,  where  p  is  the  number  of 


signals),  S  is  the  p  by  p  signal  covariance 
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where  r,  represents  the  elementary  time  delay  between  two  consecutive  sensors  for 
the  ith  arrival. 

This  formulation  allows  the  effects  of  spatial  smoothing  to  be  intei.  -eted  in  a 
useful  way  in  terms  of  the  smoothed  signal  covariance  matrix: 

S  =  77  E  S  (  *k-1  )H  (2.6) 

K  k  —  1 

where  K  is  the  number  of  averages,  also  called  the  smoothing  index.  Restricting  the 
discussion  to  two  correlated  arrivals,  [Reddy;  1987]  indicates  that  perfect 
decorrelation  is  achieved  when  the  (1,2)  term  of  S  is  equal  to  zero.  This  term  is 
given  by 

5(1.  2)  =  S  exp  \-j(k- l)2*-f  (sine „  -  sine,)]  (2.7) 

A.  k  -  1  A 


S(l.  2)  = 


sin— ——-(sinOn  -  sine,) 
gp^l  P  X  0 _ _ 

K  sin-^-(sme0-sine,) 


where  p  is  complex  correlation  between  the  two  arrivals  of  power  er^,  erf  and 
physical  arrival  angles  e0  and  e,.  d  is  the  array  interelement  spacing,  and  x  the 
acoustic  wavelength.  Perfect  decorrelation  is  achieved  when 

K  _ _ h. _ 

d  I sine0  -  sine, I 


(2.9) 


These  values  of  the  smoothing  index  K  can  be  tabulated  for  a  particular  array 
geometry  and  depends  on  the  arrival  pair  angular  separation  and  position  with 
respect  to  broadside.  For  instance,  12  averages  are  required  to  fully  decorrelate  a  5° 
separated  pair  impinging  on  an  array  with  half  wavelength  interelement  spacing  at 
10°.  Equation  (2.9)  shows  that  the  number  of  averages  necessary  to  decorrelate  two 
arrivals  depends  on  their  angular  spacing.  Two  closely  spaced  arrivals  are  very 
difficult  to  decorrelate  [Reddy;  1987]  and  require  a  large  number  of  averages  which 
in  practice  may  not  be  possible,  the  array  having  a  limited  number  of  sensors. 
Futhermore,  the  effective  aperture  is  reduced  from  M  to  s  sensor  which  results  in 
lower  resolution.  Decorrelation  is  achieved  at  the  cost  of  resolution  and 
computational  burden. 

Simulations  can  be  used  to  quantify  the  signal  cancellation  loss  at  the  output 
of  the  MV  processor.  As  before,  a  10-|-  wavelength  array  with  32  sensors  receives  a 

pair  of  coherent  unit-power  plane  wave  arrivals  in  the  set  of  incidence  angles  {  0°, 
5°,  10°,  15°,  20°,  25°,  30°,  35°,  40°}.  A  background  sensor  noise  is  also  included  with 
a  power  of  -20  dB.  The  possible  angular  spacings  are  multiples  of  5°  which  is  close 
to  the  Rayleigh  resolution  of  the  array.  For  each  pair  of  arrival  combinations,  spatial 
smoothing  followed  by  MV  beamforming  is  performed.  The  difference  in  dB 
between  the  true  power  (in  this  case  0  dB)  and  the  power  of  the  arrivals  at  the 
output  of  the  beamformer  is  measured  and  plotted  in  Figure  2.3  as  a  function  of  the 
number  of  spatial  smoothing  averages  (which  varies  from  3  to  19).  Each  panel 
corresponds  to  an  angular  separation  which  includes  a  number  of  arrival 
combinations,  e.g.  the  10°  separation  corresponds  to  the  following  pairs:  (0°,  10°). 
(5°,  15°),  (10°,  20°),  (15°,  25°).  (20°,  30°).  (25°,  35°),  (30°,  40°).  The  impact  of  the 
number  of  averages  on  the  final  angular  spectrum  can  be  observed  in  Figure  2.3  and 
a  best  subarray  length  can  be  selected  for  each  case.  Considerable  smoothing  must 
be  performed  when  the  two  arrivals  are  close  together:  a  5°  separation  requires  at 
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least  16  averages  for  the  MV  processor  to  yield  cancellation  loss  on  the  order  of  a 
few  dB.  On  the  other  hand,  a  subsegment  length  of  24  does  a  good  job  at 
decorrelating  arrivals  separated  by  at  least  10°. 

2.2.3.2.  Modified  Spatial  Smoothing 

A  second  version  of  spatial  smoothing  has  been  proposed  [Williams;  1988]  and 
uses  an  estimate  of  the  covariance  matrix  which  is  the  modified  covariance  matrix 
of  the  combined  forward/backward  linear  prediction  algorithm  for  auto-regressive 
spectral  estimation  [Burg;  1967].  The  modified  covariance  method  is  known  to 
perform  better  in  time  series  auto-regressive  spectral  estimation  in  producing  less 
spurious  peaks  [Marple;  1987].  The  modified  covariance  matrix  also  is  used  in  the 
maximum  entropy  wavenumber  processing  for  line  arrays  [Barnard;  1982]  and  in  the 
Kumaresan-Tufts  high  resolution  method  of  direction  of  arrival  estimation 
[Kumaresan;  1983]. 

In  modified  spatial  smoothing,  the  covariance  matrices  of  the  array  data  and  of 
the  reversed  data  are  estimated  and  averaged  to  produce  the  array  covariance 
matrix.  This  is  made  possible  by  the  special  geometry  of  equi-spaced  line  arrays 
and  the  planar  array  manifold.  The  modified  spatial  smoothing  does  the  average  of 
the  spatially  smoothed  covariance  matrices  estimated  on  the  99  %  overlapped 
subarrays  and  their  corresponding  reversed  matrix.  Modified  spatial  smoothing  can 
be  expressed  as  [Williams;  1988] 

M  -s 

R  smooth*  -  -M-_  \~+1  S0tt(X|  X,"  -|  J(X<X(H)*j)  (2.10) 

where  *  is  the  complex  conjugate  operation  and  J  is  the  reflection  matrix  or  reverse 
matrix  given  by 
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It  is  shown  in  [Pillai;  1989]  that  if  p  correlated  arrivals  impinge  on  the  array,  it 
is  sufficient  to  average  subarrays  of  p  sensors  to  recover  the  rank  of  the  signal 

covariance  matrix  (i.e.  to  have  det  S  ^  0),  this  result  was  first  stated  and  proved 
(under  some  specific  constraints)  in  [Williams;  1988].  With  modified  spatial 

2  M 

smoothing,  a  M  sensor  array  is  able  to  detect  at  most  correlated  arrivals 

[Williams;  1988].  This  result  is  of  great  interest  because  it  indicates  that  the  number 
of  averages,  necessary  to  determine  the  direction  of  arrivals  with  principal 
component  methods,  is  half  of  what  the  original  spatial  smoothing  requires.  It 
allows  the  so-called  increase  in  efficient  array  aperture,  thus  better  resolution. 

The  issue  with  the  MV  processor  is  the  decorrelation  rate  that  modified  spatial 
smoothing  can  achieve,  compared  to  original  spatial  smoothing.  The  study  outlined 
in  the  previous  section  is  repeated  here  for  the  modified  spatial  smoothing.  One 
considers  two  correlated  arrivals  impinging  on  a  M  sensor  array  referenced  to  its 
center  as  in  [Takao;  1987],  Then,  the  full  aperture  plane  wave  steering  vectors  are 
given  by 
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)m  m  I,  M 


(2.12) 


dJ  = 


exp (j-yy(m  -  -yy-)sin 9t 


The  array  and  signal  covariance  matrices  can  be  separated  into  a  part  that 
corresponds  to  original  spatjal  smoothing  Rlt  Si  and  a  part  that  correti,  nds  to 
modified  spatial  smoothing  R2.  S2.  The  (1,2)  term  of  the  spatially  smoothed  signal 
covariance  matrix  is  given  by  Equation  (2.8).  The  modified  part  of  the  array 
covariance  matrix  is  given  by 

R2  =  JR,*J  (2.13) 

where  Ri  the  spatially  smoothed  covariance  matrix.  Using  the  factorization  of  R,  in 
terms  of  the  signal  covariance  matrix  Si 

Ri  =  BS,Bh  (2.14) 


where  BJs  the  stacked  direction  matrix  over  the  subarrays,  the  modified  covariance 
matrix  R2  is 


R2  =  JB*Si*Br  J 


(2.15) 


Using  the  fact  that  ]B*  =  B  (Williams;  1988],  R2  =_B  S[_EH  and  since  R2  =  B  S2  BH,  the 
modified  part  of  the  signal  covariance  matrix  is  S2  =  Si'.  Therefore,  the  (1,2)  term  of 
S  is 


S(l,  2)  = 


gpgi  I/>1 
M-s+l 


sin w-^-K(sin0o  -  sin0,) 

- -j - cos  <t> 

sin jry(sm0o  -  sin0]) 


(2.16) 


where  </>  and  Ipl  are  given  by  the  complex  correlation  between  the  two  arrivals 
p  =  Ipl  eJ*. 

This  result  shows  a  peculiarity  of  modified  spatial  smoothing:  when  the 
electrical  phase  <t>  between  the  two  arrivals  referenced  at  the  center  of  the  array  is 

y,  perfect  decorrelation  is  achieved,  whatever  spatial  smoothing  length  is  used. 

The  array  covariance  matrix  is  Toeplitz  and  there  is  no  loss  due  to  correlation  at 
the  output  of  the  MV  beamformer,  even  if  s  =  M.  When  the  electrical  phase  is  zero, 
the  decorrelation  rate  is  identical  to  that  of  spatial  smoothing.  These  results,  that 
strongly  depend  on  the  electrical  phase  <f>.  are  similar  to  what  was  observed  in 
[Gabriel;  1986]  for  adaptive  antennas  and  in  [White;  1979]  for  the  minimum  entropy 
method  operating  on  fully  correlated  arrivals  in  quadrature.  The  chances  of  having 
two  signals  in  quadrature  at  the  center  of  the  array  are  slim.  Thus,  all  one  can  state 
is  that  the  modified  spatial  smoothing  works  as  well  or  better  than  the  original 
spatial  smoothing  in  decorrelating  a  pair  of  correlated  arrivals.  If  the  reference  is 
taken  at  another  spatial  position  along  the  array,  the  dependence  of  the  modified 
spatial  smoothing  performance  on  the  electrical  phase  <f>  still  exists.  The  precise 
phase  relationship  between  the  two  signals  for  which  decorrelation  is  achieved  is 

not  y  and  depends  on  the  reference  point,  the  interelement  spacing,  the 

wavelength,  and  the  arrival  angles. 

These  theoretical  results  showing  the  importance  of  the  relative  phase 
between  two  arrivals  are  simulated  with  the  same  32  sensor  array  receiving  two 
unit-power  plane  wave  arrivals  at  -5°  and  0°  incidence  angle.  The  subarray  length  is 
32  so  that  only  two  averages  using  modified  spatial  smoothing  are  performed.  The 
maximum  loss  of  the  two  arrivals  due  to  correlation  at  the  output  of  the  MV 
beamformer  is  plotted  on  Figure  2.4.  Whon  the  two  signals  are  in  quadrature,  there 
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is  no  loss  since  the  modified  covariance  matrix  is  Toeplitz.  The  loss  due  to 
correlation  increases  as  the  electrical  phase  decreases,  and  reduces  to  what  is 
produced  by  original  spatial  smoothing. 

2.2.4.  Conclusions 

The  effect  of  correlated  arrivals  on  the  adaptive  MV  beamformer  has  been 
highlighted  and  spatial  smoothing  has  been  described  and  interpreted.  The 
smoothing  techniques  attempt  to  reduce  the  modulation  of  the  diagonals  of  the 
covariance  matrix  which  should  be  Toeplitz  according  to  the  homogeneous  signal 
and  planar  wavefront  model.  Signal  cancellation  occurs  as  soon  as  correlation 
among  arrivals  exists.  The  difficulty  of  decorrelating  two  closely  spaced  arrivals 
clearly  appears  in  the  simulations.  In  practice,  substantial  loss  due  to  correlation 
may  still  occur  even  after  heavy  spatial  smoothing.  Modified  spatial  smoothing 
generally  produces  a  larger  decorrelation  than  the  original  spatial  smoothing  but 
depends  on  the  relative  phase  of  the  arrivals.  Although  spatial  smoothing 
techniques  reduce  the  effective  aperture  of  the  array  (thus  reducing  the  resolution 
and  increasing  the  computational  burden),  they  are  a  necessary  pre-processing  step 
before  using  the  adaptive  MV  processor. 

2.3.  Wavefront  Curvature 

2.3.1.  Introduction 

Besides  the  surface  and  bottom  boundaries,  another  important  aspect  of  the 
underwater  medium  is  its  highly  refractive  character,  which  is  variable  in  the 
vertical  as  well  as  the  horizontal.  Propagation  over  long  range  is  controlled  by  the 
sound  speed  profile  depth  and  range  dependence  [Murphy;  1987],  The  wavefront 
curvature  due  to  refraction  is  investigated  in  the  following. 

In  ray  theory,  the  variation  of  sound  speed  with  depth  causes  ray  bending  and 
can  be  summarized  by  Snell's  law  [Brekhovskikh;  1982]  : 

— =  constant  (2.1 7) 

where  c  is  the  local  sound  speed  and  9  the  ray  angle  with  respect  to  the  horizontal. 
Ray  bending  can  be  large,  for  instance,  a  10°  angle  ray  at  great  depth  with  speed  of 
1500  m/s  has  a  13.8°  angle  with  respect  to  the  horizontal  shallower  in  the  water 
column  where  the  sound  speed  equals  1480  m/s.  Since  ray  bending  can  be  large 
over  large  depth,  wavefront  curvature  could  have  an  important  impact  on  the 
angular  spectra  produced  by  the  beamformers  processing  the  data  from  large 
vertical  line  arrays.  Processing  a  field  with  varying  arrival  angle  across  the  array 
aperture  is  similar  to  the  case  of  processing  nonstationary  random  processes  in 
time  series  analysis.  The  selection  of  the  plane  wave  array  manifold  that 
corresponds  to  a  medium  of  a  constant  sound  speed,  may  not  be  appropriate, 
especially  for  a  very  long  vertical  line  array. 

First,  attention  is  directed  toward  generating  a  curved  wavefront  array 
manifold  based  on  the  sound  speed  profile  across  the  array.  This  manifold  will 
then  allow  the  simulation  of  a  more  realistic  pressure  field  and  the  evaluation  of 
the  plane  wavefront  beamformers  under  curved  arrivals. 

2.3.2.  Generating  Curved  Wavefronts 
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2.3.2. 1.  Derivation  of  the  Phase  Relationships  in  the  Vertical 

Wavefront  curvature  is  modelled  with  the  simple  geometric  approach  of  ray 
theory.  Since  the  sound  speed  continuously  varies  with  depth,  the  phase  variations 
from  a  point  M0  in  space  to  a  point  Mi  can  be  expressed  by  the  integral  relationship 

«1 

0  =  /F  dr  (2.18) 

M0 

where  0  is  the  phase  of  the  propagating  wavefront,  F  the  wavenumber  vector  which 
is  a  function  of  position  in  the  two  dimensional  space  sampled  by  the  vertical  array 
defined  by  the  (r,  z)  coordinate  system,  with  r=  r  cr  +  z  {srr,erz)  spans  the  space 
as  defined  by  Figure  2.5.  F  can  be  expressed  as 

F  (  r,  z)  =  ■  f  cos d(r,  z)  e;  +  sin 6{r,z)  e~z  j  (2.19) 

where  f  is  the  acoustic  frequency,  c(r  ,  z)  the  sound  speed,  9{r  ,  z)  the  ray  angle 
with  respect  to  the  horizontal  at  a  point  in  space  of  coordinates  (r  ,  z).  If  the 
receiving  array  has  its  sensors  numbered  from  top  to  bottom,  negative  angles  of 
arrival  corresponds  to  downgoing  sound  or  uplooking  beams.  Since  the  array  is 
deployed  in  the  vertical,  only  phase  variations  in  the  vertical  are  of  interest  so  that 
dr=dze%  and,  assuming  a  horizontally  stratified  ocean  (i.e.  range  independent), 
Equation  (2.18)  becomes 

0  =  2*  f  /  -f-^-dz  (2.20) 

z0  c<2> 


Using  Snell’s  law,  one  gets 


*i 

0  =  2*  f  J 

zo 


1  -  cos20(zo)  -4— r 

c(z0) 


(2.21) 


where  Sgn  is  the  signum  function,  and  cos2g(z0)[  \  ]2  <  1  is  implicitely  assumed 

for  all  z.  If  the  medium  is  isovelocity  (for  all  z,  c(z)  =  c(z0)),  Equation  (2.21) 
reduces  to  the  plane  wavefront  result : 


0  = 


2*  fsin0(zo)  . 

C(Z0— (Z'  -  Z°> 


(2.22) 


Equation  (2.21)  provides  the  variation  of  phase  in  the  vertical  and  is  identical 
to  a  ray  representation  of  the  field,  as  shown  in  [Tran;  1989a]. 


2. 3. 2. 2.  Partial  Insonification 

The  phase  relationships,  thus  the  arrival  angles,  depend  on  the  selected 
reference  depth  and  corresponding  sound  speed.  Several  reference  points  can  be 
used,  for  instance  the  center  of  the  array  aperture,  or  the  point  of  minimum  sound 
speed  across  the  array.  Using  a  particular  reference  point,  the  results  of  the 
previous  Section  can  be  used  to  compute  the  phase  entering  into  the  steering 
vectors.  The  replica  vectors  then  can  be  stacked  in  the  array  manifold  matrix  which 
has  complex  columns  that  correspond  to  physical  angles  and  rows  that  correspond 
to  the  different  sensors  of  the  array.  The  sensors  are  assumed  numbered  from  top 
to  bottom.  The  ith  sensor  or  the  ith  row  in  the  matrix  is  associated  with  a  sound 
speed  c(.  If  c,  >  c0.  where  c0  is  the  reference  sound  velocity  (sound  velocity  at  the 
center  of  the  array),  a  ray  with  angle  0O  with  respect  to  the  horizontal  at  the  center 


23 


may  not  propagate  to  the  sensor  depending  on  0O.  The  condition  for  physical 
existence  of  a  ray  at  the  ith  sensor  is 

e0  >  cos-1—  (2.25) 

Q 

If  Equation  (2.25)  does  not  hold,  the  ray  has  turned  over,  and  the  steering  vector 
complex  exponential  is  replaced  by  zero.  Beyond  the  turning  point,  there  is  a 
shadow  zone  where  the  solution  of  the  wave  equation,  under  the  WKB 
approximation,  is  exponentially  decaying  (Boyles;  1984,  p.  210].  Here,  the  ray  model 
assumes  that  there  is  no  insonification  past  the  turning  point.  Figure  2.6  shows  a 
generic  array  manifold  matrix  for  a  particular  sound  velocity  profile  with  respect  to 
a  receiving  line  array.  The  shaded  area  corresponds  to  zero-filled  entries  (i.e.  no 
sound). 

2. 3. 2. 3.  How  Curved  Are  Curved  Wavefronts  ? 

The  relative  importance  of  wavefront  curvature  with  respect  to  plane 
wavefronts  is  now  assessed.  This  is  performed  by  simply  comparing  the  phase  of  a 
curved  wavefront  to  that  of  a  corresponding  plane  wavefront  (i.e.  with  same  angle 
at  the  reference  point),  with  the  reference  taken  at  the  center  of  the  array.  A 
simulation  is  performed  on  a  64  element  array  with  25  m  spacing  in  the  water 
column  with  a  sound  speed  linearly  increasing  from  1481  m/s  at  the  top  of  the 
array  to  1515  m/s  at  the  bottom  of  the  array.  The  reference  sound  speed  at  the 
center  of  the  array  is  1498  m/s  and  the  frequency  is  20  Hz.  The  phase  differences 
for  each  sensor  position  and  all  arrivals  with  angles  between  -90°  and  90°  are 
plotted  on  Figure  2.7.  The  relative  phase  change  between  curved  and  plane 
wavefront  arrivals  is  up  to  10  %  of  the  maximum  total  phase  phase  variation. 

Figure  2.7  shows  that  only  a  few  steering  vectors  yield  a  large  phase  difference 
between  the  curved  and  plane  cases.  The  arrivals  which  are  subject  to  these  large 
deviations  have  low  angles  with  respect  to  the  horizontal.  Even  short  apertures 
experience  some  important  relative  phase  variations  at  low  angles. 

2.3.3.  Plane  Wave  Beamformers  in  a  Curved  Wavefront  Environment 

The  response  of  the  plane  wave  beamformers  to  a  curved  wavefront  field  is  of 
interest  in  order  to  assess  the  importance  of  wavefront  curvature.  The  64  element 
array  of  the  previous  Section  is  used  in  the  same  environment.  The  data  at  the 
input  of  the  plane  wave  beamformer  are  unit-power  curved  wavefront  arrivals. 
Curved  and  plane  wavefront  phase  relationships  are  generated  as  discussed  in  the 
previous  sections.  No  noise  is  included  here.  A  Kaiser-Bessel  window  with  an  a 
parameter  of  1.5  weights  the  data  across  the  array  in  order  to  yield  a  38  dB  side 
lobe  rejection.  Bearing  responses  for  arrivals  at  -40°,  -30°,  -20°,  -10°,  -5°,  -0.2°,  0.2°, 
7°,  10°,  20°,  30°.  40°  are  plotted  in  Figure  2.8.  The  bearing  responses  of  the 
conventional  beamformer  for  curved  data  are  similar  to  the  one  for  plane  wave  data 
in  the  case  of  high  angles  of  arrival.  Differences  appear  for  low  angles  of  arrival, 
where  the  main  peaks  appear  smeared  and  biased.  The  conventional  beamformer 
yields  a  peak  at  5°  with  a  3  dB  loss  under  the  0.2°  arrival.  Furthermore,  the  side  lobe 
level  for  low  arrival  angles  is  raised  by  at  least  20  dB. 

Since  there  is  no  noise  in  the  simulations,  the  MV  processor  is  extremely 
sensitive  to  mismatch  compared  to  the  conventional  processor  [Seglison;  1970, 
McDonough;  1971],  as  illustrated  in  Figure  2.9.  The  plane  wave  MV  bearing 
responses  show  large  mismatch  loss  for  all  directions  especially  for  near  horizontal 
arrivals  where  curvature  is  largest.  The  loss  due  to  mismatch  at  the  output  of  the 
MV  processor  is  plotted  in  Figure  2.10  for  different  array  lengths  in  the  same 
environment.  The  loss  increases  with  the  number  of  sensors.  It  also  increases  as  the 
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arrivals  becomes  closer  to  the  horizontal.  These  results  suggest  that  wavefront 
curvature  should  be  included  in  the  processing. 

2.3.4.  Curved  Wavefront  Beamforming 

The  results  of  the  previous  Section  indicate  that  using  the  plane  array 
manifold  for  a  vertical  array  in  a  curved  wavefront  environment  may  not  be 
appropriate,  especially  for  low  angle  of  arrivals.  Curved  wavefront  steering  vectors 
should  be  used  as  replica  vectors  [Murphy87]. 

2.3.4. 1.  Conventional  Curved  Wavefront  Beamforming 

Generally,  the  Bartlett  beamformer  is  used  with  a  tapering  window  in  order  to 
provide  good  sidelobe  rejection.  Windowing  is  necessary  to  avoid  unacceptably 
high  sidelobe  levels,  but  is  not  as  straightforward  as  in  the  plane  wave  case, 
because  of  the  partial  insonification  of  the  array  aperture.  The  window  length  has 
to  be  adjusted  and  applied  to  the  insonified  part  only  of  the  array  manifold  matrix. 
Curved  wavefront  processing  using  an  adjusted  Kaiser-Bessel  window  with  an  a 
parameter  of  1.5  avoids  the  smearing  and  bias  for  low  angle  of  arrival  but  not  the 
high  sidelobe  level  on  the  order  of  -20  dB  which  is  due  to  the  reduction  of  the 
aperture  length  corresponding  to  partial  insonification  [Tran;  1989a], 


2. 3.4. 2.  MV  Curved  Wavefront  Beamforming 


Array  Covariance  Matrix  of  a  Curved  Wavefield 

To  understand  why  the  plane  wave  MV  processor  performs  so  poorly  under 
curved  wavefront  data,  the  effects  of  curvature  on  the  array  covariance  matrix  are 
now  investigated.  Uncorrelated  plane  wavefront  arrivals  impinging  on  a  line  array 
with  equispaced  sensors  lead  to  a  Toeplitz  Hermitian  covariance  matrix.  When  the 
arrivals  have  curved  wavefronts,  the  corresponding  phase  variations  from  sensor  to 
sensor  are  not  constant  and  the  covariance  matrix  is  not  Toeplitz.  The  contribution 
to  the  covariance  matrix  of  an  arrival  with  a  signal  direction  E 


Er 


eJi  1  eJ*Z  .  .  .  e-i*M 


is  E  Eh  which  has  a  first  upper  diagonal  given  by 

fg Ml  ~  £  j^2  ...  g-^M— 1  - 


(2.26) 

(2.27) 


Because  of  curvature,  one  has  for  z> j,  0,  -  0I+1  #  <j>j  -  <f>J+l,  and  E  EH  is  not  Toeplitz. 


Spatial  Smoothing  in  a  Curved  Wavefield 

Since  correlated  arrivals  with  curved  wavefront  actually  impinge  on  the  array, 
spatial  smoothing  techniques  must  be  used.  The  direct  application  of  the  spatial 
averaging  schemes  introduces  some  distortion  since  the  underlying  structure  of  the 
covariance  matrix  is  not  Toeplitz  and  depends  on  the  arrivals  which  are  unknown  a 
priori.  Simulations  in  a  noise  free  environment  indicate  that  using  curved  wavefront 
replica  vectors  after  spatial  smoothing  still  results  in  some  loss,  although  it 
produces  an  improvement  compared  to  using  planar  wavefronts,  as  shown  in  Figure 
2.1 1  for  a  32  sensor  array  extracted  from  the  64  sensor  array  of  the  previous 
section.  This  loss  is  due  to  the  mismatch  between  the  curved  wavefront  replica 
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vectors  and  the  equivalent  signal  direction  vectors  after  spatial  smoothing.  For 
each  bearing  characterized  by  a  signal  direction  vector  E  (which  has  a  curved 
wavefront),  the  mismatch  can  be  minimized  by  choosing  the  replica  vector  E  that 
maximizes  the  power  at  the  output  of  the  conventional  processor 


Eh 


T(e  eh )  Ie 


(2.28) 


where  T  denotes  the  smoothing  transformation.  The  "best"  replica  vector  is  given 
by  the  eigenvector  of  T( E  EH)  that  has  the  largest  eigenvalue.  This  replica  vector 
will  be  referred  to  as  the  eigensteering  vector.  Because  the  distortion  effects 
associated  with  the  smoothing  transformation  are  relatively  small,  the  smoothed 
matrix  is  of  rank  1  for  all  practical  purposes. 

Figure  2.12  gives  the  bearing  responses  for  the  MV  beamformer  processing  the 
same  32  sensor  array  as  in  Figure  2.11,  after  spatial  smoothing  with  subsegment 
length  of  24.  Eigensteering  vectors  are  used  and  there  is  no  loss  due  to  mismatch. 


Performance  of  the  MV  Beamformer  after  Smoothing 

The  issue  finally  addressed  here  is  to  determine  how  the  MV  processor 
performs  under  correlated  curved  wavefront  arrivals  after  pre-processing  with  the 
spatial  smoothing  techniques.  A  simulation  (similar  to  the  plane  wave  case  of 
Section  2.2.3. 1)  with  a  pair  of  coherent  arrivals  [Tran;  1989a]  showed  that  the  MV 
processor  after  spatial  smoothing  performs  like  the  plane  wave  case  except  for  low 
angle  arrivals  near  horizontal  which  suffer  more  loss.  In  some  cases,  the  low  angle 
arrivals  are  not  even  resolved  as  a  result  of  the  partial  insonification  and  the 
reduction  of  the  efficient  aperture  because  of  spatial  smoothing  [Tran;  1989a], 

2.3.5.  Conclusions 

Based  on  the  sound  speed  profile,  curved  wavefront  steering  vectors  can  be 
derived  conveniently  using  the  simple  geometric  approach  of  ray  theory.  Wavefront 
curvature  is  large  for  near  horizontal  arrivals,  even  for  a  relatively  short  array 
aperture.  Arrivals  with  steep  angles  with  respect  to  the  horizontal  are  quite  similar 
to  plane  wave  arrivals.  These  same  results  are  observed  on  the  bearing  responses 
of  the  plane  wave  conventional  and  MV  beamformers  which  perform  reasonably 
well  for  steep  arrivals  with  respect  to  the  horizontal  but  poorly  for  near  horizontal 
arrivals.  The  great  sensitivity  to  mismatch  of  the  MV  processor  is  observed  in  the 
noise-free  simulations.  The  direct  application  of  spatial  smoothing  to  decorrelate 
correlated  curved  wavefront  arrivals  causes  some  slight  distortion  and  using  the 
curved  wavefront  array  manifold  still  yields  some  loss  due  to  mismatch.  This  loss 
due  to  mismatch  can  be  avoided  by  using  the  eigensteering  vectors.  The 
performance  of  the  curved  wavefront  MV  beamformer  under  correlated  curved 
wavefront  arrivals  is  similar  to  the  plane  wave  case  (i.e.  plane  wave  data  and  plane 
wave  replica  vectors),  except  for  near  horizontal  arrivals  where  performance 
degradation  occurs. 
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Figure  2.1:  Upper  diagonals  of  the  covariance  matrix.  Two  uncorrelated  (le 
coherent  (right)  arrivals  at  -10°  and  -20°  impinge  on  the  32  sensor  array. 
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Figure  2.4:  Signal  cancellation  loss  of  the  MV  processor  after  modified  spatial  smooth¬ 
ing.  The  loss  varies  as  a  function  of  the  electrical  phase  between  the  two  coherent  arrivals 
that  impinge  on  the  array  at  -5°  and  0°. 
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Figure  2.8:  Bearing  responses  of  the  Conventional  FFT  beamfoimer  to  curved  wave- 
fronts.  A  Kaiser-Bessel  window  a  =1.5  weights  the  data.  Panels  A  to  L  correspond  to 
arrivals  at  -40°,  -30°,  -20°,  -10°,  -5°.  -0.2°,  0.2°,  7°,  10°,  20°,  30°  and  40°. 
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Figure  2.9:  Bearing  responses  of  the  plane  wave  MV  processor  to  curved  wavefronts. 
Panels  A  to  L  correspond  to  arrivals  at  -40°,  -30°,  -20°,  -10°,  -5°,  -0.2°,  0.2°,  7°,  10°,  20°, 
30°  and  40°. 
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Figure  2.10:  Loss  due  to  mismatch  of  the  plane  wave  MV  processor  under  curved  wave- 
front  data.  M  is  the  number  of  array  sensors. 
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Figure  2.11:  Bearing  responses  of  the  curved  wavefront  MV  processor  to  curved  wave- 
fronts.  Panels  A  to  L  correspond  to  arrivals  at  -40°,  -30°,  -20°,  -10°,  -5°,  -0.2°,  0.2°,  7°, 
10°,  20°,  30°  and  40°. 
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Figure  2.12:  Bearing  responses  of  the  MV  processor  to  curved  wavefronts  using  the 
eigensteering  vectors  for  spatial  smoothing  on  24  sensors.  Panels  A  to  L  correspond  to 
arrivals  at  -40°,  -30°,  -20°,  -10°,  -5°,  -0.2°,  0.2°,  7°,  10°,  20°,  30°  and  40°. 
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3.  Processing  Synthetic  Data 


3.1.  Overview 

Full  and  subaperture  beamforming  approaches  are  proposed  to  process  the 
data  from  large  vertical  line  arrays  and  for  the  purpose  of  estimating  signal  arrival 
structure.  Simulations  on  synthetic  data  produced  by  a  normal  mode  model  are 
used  to  determine  if  wavefront  curvature  must  be  considered  and  if  the  processing 
approaches  appear  qualitatively  and  quantitatively  viable.  Simulations  involving  a 
3000  m  long  vertical  line  array  indicate  that  curved  wavefront  replica  vectors  based 
on  the  sound  speed  profile  should  be  used  to  estimate  the  vertical  angular 
spectrum,  especially  if  the  arrival  structure  contains  low  arrival  angles.  Then, 
simulations  involving  a  900  m  aperture  and  128  sensors  show  that  subaperture 
processing  is  workable  and  potentially  worthwile.  For  such  an  aperture,  it  is  shown 
that  plane  wavefront  replica  vectors  can  be  used  without  any  concern  of  mismatch, 
even  at  high  input  signal-to-noise  ratio,  the  signal-to-noise  ratio  at  a  single  element 
(i.e.  10  dB).  Processing  gain  and  resolution  performance  relative  to  the  conventional 
processor  are  quantified  for  various  input  signal-to-noise  ratio  scenarios.  Using  the 
Minimum  Variance  processor  results  in  some  reduction  of  the  processing  gain.  This 
reduction  decreases  with  lower  input  signal-to-noise  ratios,  as  the  sensitivity  to 
mismatch  of  the  adaptive  processor  decreases.  Resolution  enhancements  can  be 
obtained  with  the  Minimum  Variance  processor  at  high  and  medium  input  signal-to- 
noise  ratios  (i.e.  10,  -10  dB),  but  not  for  lower  input  signal-to-noise  ratios. 

3.2.  Introduction 

This  research  investigates  how  to  process  the  data  from  large  vertical  line 
arrays  such  as  the  Marine  Physical  Laboratory  (MPL)  digital  acoustic  array  [Sotirin; 
19891,  a  120  sensor  array  cut  for  100  Hz  and  sampling  900  m  of  water  column.  As 
indicated  in  the  previous  chapters,  the  focus  is  put  on  how  to  use  the  adaptive 
Minimum  Variance  processor  in  practice.  The  MV  processor  can  provide  substantial 
resolution  improvement  but  lacks  robustness  and  exhibits  a  large  sensitivity  to 
mismatch  at  high  input  signal-to-noise  ratios. 

Full  aperture  coherent  processing  offers  the  best  processing  gain  in  theory.  It 
is  commonly  performed  using  the  conventional  FFT  beamformer,  which  is  a  very 
computationally  efficient  processor  [De  Fatta;  1988].  The  FFT  beamformer  does 
across-sensor  FFTs  on  a  selected  narrowband  frequency  to  obtain  angular  spectra 
that  are  then  averaged  to  achieve  statistical  stability.  Full  aperture  processing  with 
the  Minimum  Variance  beamformer  on  a  long  and  well-filled  array,  (i.e.  with  the 
number  of  sensors  on  the  order  of  a  hundred),  can  be  performed,  although  several 
aspects  do  not  appear  practical.  As  a  preliminary  step,  the  adaptive  processor 
requires  the  estimation  of  the  array  covariance  matrix  and  needs  a  large  number  of 
time  snapshots  in  order  to  correctly  estimate  an  invertible  covariance  matrix  and  to 
control  the  power  spectrum  variance.  The  amount  of  time  snapshots  required  by 
the  MV  processor  may  be  well  beyond  the  length  of  the  data  stationary  sample 
record,  especially  if  the  signals  of  interest  come  from  a  moving  source. 
Furthermore,  the  MV  processor  involves  matrix  operations  for  which  computation 
counts  have  a  N3  order  where  N  is  the  size  of  the  matrix  [Rung;  1988).  Large 
computational  ressources  may  be  required. 

Another  strategy  is  to  process  subarrays  with  the  MV  processor  and 
incoherently  recombine  the  corresponding  angular  spectra,  as  it  is  commonly  done 
in  time  series  spectral  estimation.  This  approach  presents  several  advantages:  it 
requires  less  snapshots  to  estimate  the  smaller  size  subarray  covariance  matrices 
and  less  computing  ressources.  Although  the  MV  processor  operates  o  i  smaller 
apertures,  it  could  yield  sufficient  resolution,  if  not  better  resolution  than 
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conventional  beamforming,  depending  on  the  input  signal-to-noise  ratio  and  should 
enjoy  a  certain  level  of  robustness  as  a  result  of  incoherent  averaging.  In  addition, 
this  approach  produces  angular  spectra  corresponding  to  different  windows  in  the 
water  column  and  may  provide  some  information  on  multipath  variability  in  the 
vertical,  if  any. 

The  implementation  of  the  two  envisioned  processing  approaches  is  discussed 
below  and  the  techniques  are  validated  using  synthetic  data.  The  simulated  data  are 
created  with  an  acoustic  normal  mode  model  that  outputs  the  complex  wavefield  at 
a  particular  frequency  for  a  particular  source-receiver  geometry.  Normal  mode 
models  provide  a  complete  solution  of  the  wave  equation  including  full  diffraction 
effects  [Boyles;  1984,  p.  199]  and  should  produce  a  realistic  wavefield  to  process 
with  the  proposed  methods.  A  deep  water  CW  tone  propagation  experiment  is 
simulated  with  two  receiving  array  geometries.  The  first  one  is  of  a  very  long 
vertical  line  array  with  an  aperture  length  over  3000  m.  The  second  geometry 
involves  a  large  vertical  line  array,  with  a  length  on  the  order  of  900  m,  comparable 
to  the  MPL  digital  array. 

3.3.  Beamformer  Implementation 

The  conventional  beamformer  and  the  adaptive  MV  processor  described  in 
Chapter  1,  are  used  in  the  following.  The  conventional  beamformer  operates  on  the 
full  aperture.  The  corresponding  angular  spectrum  will  be  used  as  a  reference.  The 
MV  processor  processes  the  full  aperture  and  subapertures  with  various  spatial 
smoothing  lengths.  Two  wavefront  models  are  considered,  the  plane  wavefront  and 
the  curved  wavefront  based  on  the  sound  speed  profile  across  the  array,  as 
developed  in  Chapter  2.  The  angular  spectra,  produced  by  the  various  processors 
using  the  two  wavefront  models,  are  referenced  to  a  particular  depth  and  sound 
speed  in  order  to  be  intercompared.  Here,  the  reference  depth  is  taken  at  the 
sound  speed  minimum  across  the  whole  water  column  (i.e.  the  sound  axis  depth). 
The  receiving  arrays  are  deployed  across  the  sound  axis  at  a  known  absolute  depth. 

3.3.1.  Using  Planar  Wavefonts 

The  planar  wavefronts,  referenced  to  the  sound  axis  depth  and  sound  speed, 
are  generated  using  Equation  (2.22).  The  differential  phase  across  the  array 
depends  only  on  the  reference  sound  speed  and  the  reference  depth  just  introduces 
an  unimportant  phase  shift.  This  is  consistent  with  the  FFT  implementation  of  the 
conventional  processor  (c.f.  Equation  (1.12))  which  produces  a  power  spectrum  in 
normalized  wavenumber  or  electrical  angle,  then  transformed  into  physical  angles 
through  a  nonlinear  mapping  (given  by  Equation  (1.9))  involving  only  the  reference 
sound  speed.  A  Kaiser-Bessel  window  with  a  parameter  of  1.5  is  used  to  weight  the 
aperture,  prior  to  the  conventional  beamforming  to  control  the  sidelobe  level  to  -35 
dB.  The  MV  processor  is  directly  implemented  with  the  quadratic  form  in  Equation 
(1.36).  In  the  case  of  processing  the  full  aperture  with  planar  wavefronts,  the  replica 
vectors  are  created  according  to  the  efficient  aperture  length  after  spatial 
smoothing  using  the  sound  axis  sound  speed  as  reference. 

Subaperture  processing  requires  the  computation  of  a  reference  sound  speed 
for  each  subaperture  in  order  to  compute  its  associated  plane  wave  replica  vectors. 
The  reference  sound  speed  can  be  the  mid-subarray  sound  speed,  the  arithmetic  or 
harmonic  sound  speed  average  across  the  subarray.  Because  the  relative  difference 
in  sound  speed  is  small  anyway,  any  choice  appears  acceptable  for  producing  the 
planar  replica  vectors.  The  harmonic  sound  speed  is  used  here  [Polvani;  1984].  Since 
different  averaged  sound  speeds  are  associated  with  each  subaperture,  a  first  order 
correction  is  necessary  to  take  into  account  ray  bending  and  insure  that  the  arrivals 
are  properly  aligned  from  one  subaperture  to  another.  The  mapping  is  performed 
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by  applying  a  correction  based  on  Snell’s  law  which  shifts  and  aligns  the  peaks  of 
the  angular  spectra  which  then  can  be  incoherently  averaged  to  produce  a 
composite  angular  spectrum.  For  example,  assume  a  two  non-overlapping 
subaperture  case  where  the  upper  subaperture  has  a  lower  reference  sound  speed 
than  the  the  lower  subaperture.  Then,  denote  Pup(0up)  and  Piow{Olow)  the  angular 
spectra  of  the  upper  and  lower  apertures.  Ptow  is  remapped  for  a  sound  speed 
deviation  from  qow  to  cup  using  Snell's  law 

cos (8,ow)  _  cos(eup) 

QOW  Oyp 

Since  Ptow  is  evaluated  on  a  regular  grid  in  0tow  angles,  which  actually  correspond  to 
an  irregular  grid  in  9up  angles  referenced  to  cup,  one  just  needs  to  do  an 
interpolation  to  go  from  an  irregular  angular  grid  back  to  a  regular  grid.  This 
procedure  is  summarized  in  Figure  3.1.  The  result  for  Piow  is  large  angular  shifts  for 
low  angles  and  small  angular  shifts  for  steep  angles  with  respect  to  the  horizontal. 
The  low  angle  arrivals  for  the  lower  aperture  corresponds  to  steeper  angle  arrivals 
for  the  remapped  lower  aperture  and  very  low  angle  sound  may  not  exist.  One  must 
note  at  this  point  that  small  variations  in  the  subarray  reference  sound  speeds  may 
result  in  large  shifts  in  angles,  especially  for  low  arrival  angles  (below  10°).  The 
harmonic  sound  speed  was  found  to  produce  a  very  good  alignment  of  the 
subaperture  peaks. 


3.3.2.  Using  Curved  Wavefronts 

Curved  wavefront  beamforming,  proposed  in  Chapter  2,  is  summarized  here. 
The  sound  speed  profile,  assumed  known  across  the  array,  is  used  in  Equation 
(2.21)  to  create  the  curved  wavefront  replica  vectors  corresponding  to  the  full 
aperture.  If  the  ray  corresponding  to  a  particular  angle  at  the  sound  axis  has  turned 
over  and  does  not  insonify  part  of  the  array,  the  non-insonified  sensors  or  rows  in 
the  array  manifold  matrix  have  zero  entries.  When  the  Bartlett  processor  beamforms 
the  data,  a  Kaiser-Bessel  window  with  a  parameter  of  1.5  is  used  to  weight  the 
insonified  part  of  the  array  like  in  the  plane  wave  case.  Care  is  necessary  to  ensure 
proper  normalization  of  the  power  spectrum,  if  the  top  sensors  and  the  bottom 
n2  sensors  are  not  insonified  by  a  particular  steering  vector  Ec,  the  curved 
wavefront  conventional  processor  is 

Pmser  =  (w  Ec)hR  (W  Ec)  (3.2) 

where  Ec  and  R  are  defined  as  before,  W  is  defined  by 
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where  the  tf,’s  are  the  Kaiser-Bessel  window  coefficients  weighting  the  M,  insonified 
sensors  of  the  array,  =  M  -  nx  -  n2  . 

When  the  MV  processor  is  used  to  process  the  full  aperture,  the  eigensteering 
vector  corresponding  to  each  curved  wavefront  replica  vector  must  be  computed  to 
avoid  mismatch  due  to  spatial  smoothing.  The  array  covariance  matrix 
corresponding  to  each  replica  vector  is  passed  through  spatial  smoothing,  the 
smoothed  matrix  is  decomposed  in  eigenvalues  and  eigenvectors  and  the  unit-norm 
eigenvector  Ee  corresponding  to  the  largest  eigenvalue  is  stored  in  the  eigen  array 
manifold.  Therefore,  a  large  number  of  computations  is  required  for  each 
smoothing  scenario.  As  in  the  case  of  the  Bartlett  processor,  care  is  necessary  to 
correctly  normalize  the  power  spectrum.  The  Minimum  Variance  power  spectrum  is 
given  by 


Pmv 


1 _ 1 

*i(E«)  EeH  T(R)’1  Ee 


(3.5) 


where  T(R)S  is  the  spatially  smoothed  covariance  matrix  on  s-sensor  subarrays  and 
S/(Ee)  is  the  number  of  sensor  actually  insonified  by  an  arrival  with  replica  vector 
Ee- 

When  subapertures  are  processed,  the  part  of  the  full  aperture  curved 
wavefront  replica  vectors  that  corresponds  to  a  particular  subaperture  is  used  to 
create  the  corresponding  eigensteering  vectors.  The  subaperture  angular  spectra 
then  are  computed  as  in  the  full  aperture  MV  angular  spectrum  using  their 
respective  eigen  array  manifold.  The  subaperture  spectra  can  then  be  directly 
averaged  together  since  they  correspond  to  the  same  depth  and  sound  speed 
references.  Therefore,  there  is  an  interesting  benefit  of  using  curved  wavefront 
replica  vectors  compared  to  plane  wavefront  replica  vector.  The  remapping 
operation  of  the  angular  spectra  (associated  with  a  particular  subarray  and  a 
particular  constant  sound  speed),  which  inherently  is  sensitive  to  the  subaperture 
sound  speeds,  is  avoided.  The  disadvantage  is  the  large  computational  complexity 
since  an  eigenvalue-eigenvector  decomposition  is  required  for  each  beam  in  order 
to  compute  the  eigen  array  manifold. 


3.4.  Processing  a  Very  Large  Vertical  Line  Array 


3.4.1.  Modeling  the  Environment 

The  simulated  data  created  by  the  ATLAS  normal  mode  model  [Gordon;  1985], 
corresponds  to  a  range  independant  deep  water  environment  with  a  water  depth  of 
4667  m.  The  sound  speed  profile,  plotted  in  Figure  3.2,  is  characteristic  of  the 
North-East  Pacific,  with  a  sound  axis  around  500  m.  The  critical  depth  of  this 
oceanic  waveguide  is  on  the  order  of  3313  m.  Thus,  the  sound  channel  has  a  depth 
excess  and  convergence  zone  propagation  may  take  place.  The  bottom  is 
characterized  by  a  frequency  invariant  bottom  loss  table  (bottom  loss  in  dB  as  a 
function  of  incidence  angle),  given  in  Table  3.1,  which  corresponds  to  a  Generic 
Sonar  Model  (GSM)  Province  Type  3  [Weinberg;  1985].  A  5  m  deep  omnidirectional 
source  is  projecting  a  20  Hz  CW  tone  with  a  source  level  of  200  dB  re  1  /iPa  at  1  m. 
The  environment  is  assumed  to  be  noise  free.  The  very  long  vertical  line  array  is 
cut  for  30  Hz  with  an  interelement  spacing  of  25  m,  has  128  equi-spaced  sensors 
and  extends  across  most  of  the  water  column  from  100  m  to  3275  m.  It  has  its 
sensors  numbered  from  top  to  bottom  so  that  at  the  output  of  the  beamformers 
positive  angles  of  arrival  correspond  to  downlooking  beams  or  upgoing  sound. 

This  receiving  array  is  assumed  deployed  56  km  from  the  source  in  its  first 
convergence  zone.  The  geometry  is  summarized  in  Figure  3.3.  The  analysis  and 
interpretation  of  these  synthetic  data  [Tran;  1989b]  are  performed  with  the  ATLAS 
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norma!  mode  mode!  and  the  GSM  ray  model  [Weinberg;  1985]  using  the  continuous 
gradient  approximation  of  the  sound  speed  profile  (CONGRATS).  The  GSM  is  used 
to  cross-check  the  ATLAS  results.  A  fair  agreement  was  found  between  the  two 
models,  although  ihe  ray  model  reports  a  well  defined  convergence  zone  around  56 
km,  while  the  normal  mode  model  does  not.  The  GSM  produces  ray  diagrams  which 
can  be  compared  to  the  beamformer  outputs.  The  eigenrays  connecting  the  source 
to  the  array  receivers,  1  ’tween  -45°  and  45°,  are  plotted  in  Figure  3.4.  The  sign 
convention  of  the  GSM  eigenrays  is  the  same  as  for  the  beamformers,  positive  angle 
eigenrays  are  propagating  upward  (i.e.  the  signs  are  flipped  with  respect  to  the 
original  Generic  Sonar  Model  results).  The  effect  of  ray  bending  with  depth  is 
clearly  observed  for  the  rays  with  smaller  angles  with  respect  to  the  horizontal. 


3.4.2.  Processing  Results 


3.4.2. 1.  Conventional  Beamforming 

The  angular  spectra  produced  by  the  conventional  beamformer  using  the  plane 
and  curved  wavefront  replica  vectors  are  overlaid  in  Figure  3.5  with  the  GSM 
eigenrays  at  the  sound  axis.  Since  no  noise  is  included  in  the  simulation,  over 
twelve  multipath  arrivals  are  resolved,  mostly  between  -30°  and  30°.  Using  either 
curved  or  planar  wavefront  replica  vectors  gives  similar  results  but  some  arrivals 
are  reported  at  slightly  different  angles,  even  with  the  reference  taken  at  the  sound 
axis  in  all  cases.  The  GSM  eigenray  structure  at  the  sound  axis  is  similar  to  the 
angular  spectra  although  the  eigenray  peak  powers  are  higher  than  the  ATLAS 
spectrum  peaks.  This  is  consistent  with  the  modeling  results.  When  curved 
wavefront  replica  vectors  are  used,  the  two  arrivals  near  -15°  and  20°  are  better 
resolved  and  their  directions  of  arrival  are  closer  to  the  GSM  eigenray  angles.  The 
conventional  beamformer  using  curved  wavefront  replica  vectors  reports  high 
power  levels  at  near  horizontal  angles.  This  is  due  to  the  partial  insonification  of 
the  aperture  which  results  in  higher  sidelobes  and  lower  resolution.  The  insonified 
part  of  the  128  element  array  is  plotted  in  Figure  3.6  as  a  function  of  arrival  angle. 
Near  horizontal  arrivals  insonify  a  very  limited  part  of  the  aperture,  if  not  a  single 
sensor  (#  13).  Since  the  beamwidth  of  a  few  sensor  array  steered  in  a  near 
horizontal  direction  is  large,  power  corresponding  to  steeper  arrivals  is  reported  at 
low  angle  of  arrival. 


3. 4. 2. 2.  MV  Processing 

This  simulation  does  not  include  noise  and  therefore  corresponds  to  an 
"infinite"  signal-to-noise  ratio.  In  such  conditions,  the  MV  processor  exhibits  a  large 
sensitivity  to  mismatch,  and  only  curved  wavefront  replica  vectors  are  used  in  the 
following.  The  angular  spectrum  produced  by  the  full  aperture  MV  beamformer 
after  spatial  smoothing  on  64  sensors  is  overlaid  in  Figure  3.7  with  the  GSM 
eigenrays  and  the  conventional  processor  using  curved  wavefront  replica  vectors.  A 
large  amount  of  smoothing  is  necessary  to  decorrelate  a  cluster  of  over  a  dozen 
multipaths.  The  MV  peak  powers  of  the  arrivals  are  up  to  20  dB  lower  than  the  ones 
produced  by  the  conventional  processor.  Since  a  large  amount  of  spatial  emoothing 
has  been  performed  and  eigensteering  vectors  are  utilized,  this  loss  is  thought  to 
be  the  result  of  an  insufficient  sampling  in  angle  (only  256  beams  are  used  between 
-90°  and  90°).  Improvement  in  resolution  is  observed  clearly  compared  to  the 
conventional  processor,  even  with  the  efficient  smoothed  aperture  reduced  by  half. 
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As  in  the  case  of  the  conventional  curved  wavefront  processor,  an  increase  in  power 
near  horizontal  angles  is  due  to  the  partial  insonification  of  the  equivalent 
smoothed  aperture,  as  shown  in  Figure  3.8  where  the  insonified  efficient  aperture  is 
plotted. 

The  processing  results,  produced  by  the  MV  processor  on  32-sensor  subarrays, 
after  spatial  smoothing  on  20  sensors  and  with  the  eigensteering  vectors,  are 
presented  in  a  waterfall  plot  in  Figures  3.9.  The  angular  spectra  correspond  to  48 
subarrays  sliding  from  the  top  to  the  bottom  of  the  full  aperture  by  two  elements 
(i.e.  the  first  subarray  extends  from  sensor  1  to  32,  the  second  subarray  extends 
from  sensor  3  to  34  and  so  on...)  and  are  all  referenced  to  the  sound  axis.  The  flat 
area  at  low  angles  corresponds  to  an  angular  region  where  no  arrival  or  sound  can 
physically  exists,  according  to  the  ray  geometric  approach  selected  here.  It  is  the 
result  of  the  partial  insonification  discussed  earlier.  The  peak  power  estimates  are 
lower  than  the  ones  reported  by  the  conventional  processor  because  of  signal 
cancellation  (the  subarrays  are  shorter  and  allow  only  limited  multipath 
decorrelation),  and  unsufficient  angular  sampling  (as  in  the  case  of  full  aperture 
processing,  only  256  beams  are  used).  A  great  variability  in  the  detected  arrivals  is 
observed  from  subarray  to  subarray,  even  though  all  the  spectra  are  referenced  to 
the  sound  axis.  It  is  largest  in  the  upper  half  and  the  middle  of  the  128  sensor 
aperture.  As  a  result,  the  incoherent  averaging  of  the  subarray  angular  spectra 
smears  the  peaks  in  the  spectrum  [Tran;  1989b]. 


3.4.3.  Conclusions 

These  simulations  on  a  very  large  vertical  line  array  indicates  that  using  curved 
wavefront  replica  vectors  with  the  conventional  processor  can  enhance  the 
resolution  of  some  arrivals,  especially  the  ones  with  small  angles  with  respect  to 
the  horizontal.  Full  aperture  MV  processing  using  the  eigensteering  vectors 
produces,  as. expected,  a  resolution  enhancement  in  this  pure  signal  environment. 
Subaperture  processing  with  the  MV  processor  and  the  eigensteering  vectors 
appears  to  be  viable  scheme,  it  indicates  the  large  variability  with  depth  of  the 
multipath  arrivals  which  characterizes  the  particular  simulated  geometry.  Thus, 
averaging  incoherently  the  subapertures  does  not  necessarily  enhance  the  peaks. 
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3.5.  Processing  A  Large  Vertical  Line  Array 


3.5.1.  Modeling  Results 

The  simulated  environment  is  the  same  as  in  the  previous  Section.  But  the 
propagation  experiment  now  involves  a  shorter  aperture  array  with  128  sensors,  cut 
for  107  Hz  with  an  interelement  spacing  of  7  m.  The  vertical  line  array  is  deployed 
from  300  to  1189  m.  The  source  is  deployed  at  5  m  depth  and  56  km  range,  and 
projects  a  100  Hz  CW  tone  with  a  level  of  200  dB  re  /xPa  at  1  m.  The  geometry  of 
the  propagation  experiment  is  summarized  in  Figure  3.3.  As  before,  the  acoustic 
modeling  is  performed  with  the  ATLAS  normal  mode  model  and  the  GSM  ray  model 
[Tran;  1989b].  Both  models  indicate  a  convergence  zone  between  50  and  60  km, 
which  is  better  defined  in  the  ray  solution  than  in  the  normal  mode  results.  The 
GSM  eigenrays  that  connect  the  source  to  the  array  receivers  are  plotted  in  Figure 
3.10,  where  positive  angle  eigenrays  are  upgoing.  The  eigenray  arrival  structure 
does  not  vary  significantly  over  depth  as  in  the  case  of  the  very  large  array  which 
was  three  times  larger. 


3.5.2.  Processing  Results 


3.5.2. 1.  Conventional  Beamforming 

The  array  aperture  is  smaller,  hence  the  effects  of  the  wavefront  curvature 
have  a  lesser  impact.  The  angular  spectra  referenced  to  the  sound  axis  produced  by 
the  conventional  beamformer,  using  planar  and  curved  wavefront  replica  vectors 
are  overlaid  in  Figure  3.11  with  the  GSM  eigenrays.  The  angular  spectra  are 
identical  whether  curved  or  plane  wavefront  replica  vectors  are  used,  except  for  the 
spurious  horizontal  arrival  due  to  the  partial  insonification  of  the  aperture  by 
curved  wavefront  arrivals.  The  0°  curved  arrival  insonifies  only  three  sensors  [Tran; 
1989b].  The  ATLAS  angular  spectra  closely  match  the  GSM  eigenray  structure.  In 
addition  to  the  good  fit  of  the  directions  of  arrival,  the  absolute  power  levels  of  the 
ATLAS  spectral  peaks  and  the  GSM  eigenrays  are  similar. 


3. 5. 2. 2.  MV  Processing 

The  angular  spectrum  produced  by  the  MV  beamformer  processing  the  full 
aperture  and  using  eigensteering  vector  that  correspond  to  a  spatial  smoothing  on 
64  sensors,  is  overlaid  in  Figure  3.12  with  the  spectrum  of  the  conventional 
processor  using  curved  wavefront  replica  vectors  and  the  GSM  eigenrays.  In  this 
pure  signal  simulation,  using  the  adaptive  processor  yields  much  better  resolution. 
The  peak  power  levels  produced  by  the  adaptive  processor  are  up  to  20  dB  lower 
than  the  conventional  processor  levels.  As  mentioned  already,  this  is  probably 
mismatch  loss  resulting  from  an  insufficient  0.7°  sampling  in  angle.  The  slight 
increase  of  power  at  low  arrival  angles  corresponds  to  the  partial  insonification  of 
the  equivalent  smoothed  aperture,  only  32  sensors  of  the  equivalent  smoothed 
aperture  are  insonified. 

The  results  of  processing  32-sensor  subapei  tures  with  the  MV  method  using 
the  eigensteering  vector  for  a  spatial  smoothing  on  20  sensors,  are  presented  in 
Figure  3.13.  These  48  spectra,  corresponding  to  32  sensor  subarrays  sliding  by  two 
sensors  from  the  top  to  the  bottom  of  the  full  aperture,  are  referenced  to  the  sound 
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axis.  The  reduction  of  the  aperture  limits  the  amount  of  spatial  smoothing,  hence, 
the  number  of  peaks  in  Figure  3.13  is  roughly  half  of  what  is  detected  with  full 
aperture  coherent  processing.  Mismatch  is  apparent  with  peak  levels  up  to  20  dB 
below  the  conventional  processor  peak  levels.  The  stronger  and  lower  angle  arrivals 
are  stable  with  depth  and  they  are  enhanced  when  the  incoherent  average  of  the  48 
angular  spectra  is  computed.  The  composite  angular  spectrum  is  plotted  in  Figure 
3.14  and  has  low  angle  arrival  'evels  on  the  same  order  as  the  full  aperture  MV 
processor  levels. 


3. 5. 2. 3.  Conclusions 

These  simulations  show  that  for  a  shorter  aperture  and  for  a  typical  arrival 
structure  with  arrivals  steeper  than  10°,  there  is  no  need  to  use  curved  wavefronts 
with  the  Bartlett  (conventional)  processor.  Beamforming  the  full  aperture  with  the 
adaptive  processor  and  curved  wavefronts  produces  a  large  increase  in  resolution 
compared  to  conventional  beamforming.  Subaperture  processing  seems  to  perform 
as  indicated  in  Section  3.2.  The  main  features  of  the  arrival  structure  are  stable  with 
depth,  therefore  incoherent  averaging  produces  some  processing  gain. 


3.5.3.  Signal-Plus-Noise  To  Noise  Ratio  Performance  Study 


3.5.3. 1.  Introduction 

Until  now,  no  component  of  noise  has  been  included.  Since  the  MV  processor  is 
very  sensitive  to  mismatch  in  a  high  signal-to-noise  ratio  environment,  only  curved 
wavefronts  (modified  to  account  for  spatial  smoothing)  have  been  used.  The 
conventional  beamforming  results  are  identical  whether  or  not  plane  or  curved 
wavefront  replica  vectors  are  used.  Thus,  one  can  question  if  using  curved 
wavefront  replica  vectors  is  actually  necessary.  A  way  to  answer  this  question  is  to 
measure,  for  a  realistic  input  signal-to-noise  ratio  (SNR)  the  mismatch  loss  (if  any) 
when  curvature  is  neglected.  This  is  performed  through  a  peak  level  performance 
study,  which  will  also  indicate  the  relative  processing  gain  performance  of  the 
various  processors.  This  study  is  performed  only  on  the  900  m  aperture  case  in 
preparation  to  processing  real  data  from  the  MPL  digital  array  (Chapter  4). 


3. 5. 3. 2.  Description  of  the  Simulation 

The  synthetic  pressure  field,  created  by  the  ATLAS  model,  is  embedded  in  a 
spatially  white  noise  or  sensor  noise  with  three  different  input  signal-to-noise  ratio 
(SNR):  -18  dB,  -1 1  dB  and  9  dB.  A  spatially  white  noise  structure  corresponds  to  a 
noise  field  dominated  by  strong  surface  generated  noise  produced  by  high  wind 
speed  conditions.  To  measure  the  potential  mismatch  caused  by  using  planar 
wavefront  replica  vectors,  the  number  of  beams  in  this  study  is  increased  to  yield  a 
very  fine  sampling  in  angle.  A  0.022°  quantization,  which  corresponds  to  4096 
beams  between  -45°  and  45°  was  found  sufficient  [Tran;  1989c].  For  each  input  SNR, 
the  field  is  processed  by  the  conventional  FFT  beamformer,  the  MV  processor 
operating  on  the  full  aperture  after  spatial  smoothing  on  92,  78  or  64  sensors,  and 
finally  by  the  MV  processor  operating  on  two  non-overlapping  subpertures  after 
spatial  smoothing  on  48  or  32  sensors.  For  each  processing  scheme,  the  peak  level 
or  output  signal-plus-noise  level  of  each  arrival  is  measured  relative  to  the  output 
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noise  level.  This  relative  peak  level  is  a  biased  measure  of  the  output  SNR  at  low 
input  SNR,  but  is  a  reliable  indicator  of  the  output  SNR  at  high  input  SNR.  Since  in 
this  simulation  the  covariance  matrices  are  perfectly  known,  very  low  level  arrivals 
with  negative  output  signal-to-noise  ratio  can  be  detected  against  the  flat 
background  noise,  and  are  included  in  the  analysis  for  completeness,  although  in  a 
real  situation,  these  arrivals  would  be  covered  by  the  noise  fluctuations.  The 
signal-plus-noise  to  noise  ratios  of  the  arrivals  detected  with  the  various  methods 
are  compared  to  the  corresponding  arrival  peak  levels  produced  by  the 
conventional  FFT  beamformer.  The  results  are  summarized  in  Figures  3.15,  3.16  and 
3.17  for  the  -18,  -1 1  and  9  dB  input  SNR. 


3. 5. 3. 3.  Results 


Plane  Versus  Curved  Wavefronts 

A  first  inspection  of  the  simulation  results  shows  that  there  is  no  noticeable 
difference  whether  planar  or  curved  (modified  to  account  for  spatial  smoothing) 
wavefront  replica  vectors  are  used.  There  is  no  mismatch  associated  with  planar 
wavefront  replica  vectors,  even  for  high  input  SNR,  for  this  typical  arrival  structure 
and  900  m  aperture.  Thus,  plane  wavefronts  can  be  used  with  the  Minimum 
Variance  processor. 


Full  Aperture  Processing  Analysis 

The  output  signal-plus-noise  to  noise  ratio  (relative  peak  level)  produced  by 
the  MV  processor  is  generally  below  the  conventional  processor  relative  peak  level. 
This  loss  increases  with  increasing  relative  peak  levels  which  is  not  surprising  since 
the  sensitivity  of  the  MV  processor  to  mismatch  increases  with  SNR.  At  low  and 
medium  input  SNR,  the  reduction  in  output  signal-plus-noise  to  noise  ratio  is  on  the 
order  of  half  a  dB  and  a  dB,  respectively.  The  multipaths  have  been  well 
decorrelated  by  spatial  smoothing,  since  in  those  two  cases,  there  is  more  loss 
when  the  spatial  smoothing  length  is  decreased.  The  high  input  SNR  case  has  the 
same  behaviour,  except  for  arrivals  above  20  dB,  where  the  relative  peak  level 
reduction  is  largest,  on  the  order  of  8  dB,  for  a  smoothing  length  of  92.  That  loss  is 
reduced  to  3  dB  when  larger  smoothing  lengths  are  used,  which  indicates  a  signal 
cancellation  limited  situation.  The  largest  arrival  around  15°  actually  corresponds 
to  a  cluster  of  arrivals  (c.f.  Figure  3.11). 


Subaperture  Processing 

Processing  two  subapertures  with  the  MV  beamformer  produces  a  reduction  of 
the  output  signal-plus-noise  to  noise  ratio,  compared  to  the  Bartlett  processor 
operating  on  the  full  aperture.  The  loss  is  on  the  order  of  half  a  dB  for  low  level 
arrivals  (i.e.  arrivals  with  a  relative  peak  level  below  2  dB)  and  on  the  order  of  a  few 
dB  for  medium  level  arrivals  (arrivals  with  relative  peak  levels  between  3  and  10 
dB),  and  up  to  6  dB  for  large  level  arrivals  (arrivals  with  relative  peak  levels  above 
10  dB). 
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3. 5. 3.4.  Array  Gain  Calculations 

Assuming  that  the  detection  of  the  signals  is  based  on  the  highest  peak  level 
arrival,  one  can  estimate  the  array  gain  achieved  by  the  array  processed  by  the 
various  beamformers.  The  array  gain  is  given  by  the  signal-to-noise  ratio  at  the 
output  of  the  beamformer  divided  by  the  input  signal-to-noise  ratio  at  the 
individual  phones.  The  input  signal-to-noise  ratios  are  fixed  in  the  simulation.  The 
highest  output  signal-to-noise  ratio  can  be  easily  derived  from  the  highest  relative 
peak  level  measured  in  the  angular  spectra  by  substracting  the  contribution  of  the 
noise.  The  estimates  of  the  array  gain  are  plotted  in  Figure  3.18  for  the  various 
methods.  They  are  always  lower  than  10  log  M,  or  21  dB,  the  theoretical  array  gain 
assuming  uncorrelated  arrivals  of  a  known  wavefront  in  an  incoherent  noise  field 
[Urick;  1983,  p.  37,  Ziomek;  1985,  p.  121].  The  array  gain  of  the  Bartlett  processor  is 
constant  with  input  SNR  and  equal  to  14  dB.  For  low  and  medium  SNRs  (-18  and  -11 
dB),  the  array  gain  of  the  MV  processor  is  close  to  the  Bartlett  processor  array  gain, 
and  essentially  the  same  for  the  low  and  medium  input  SNRs  (-18  and  -11  dB),  with 
a  reduction  on  the  order  of  2  dB  and  3.5  dB  when  processing  two  subarrays  with 
spatial  smoothing  lengths  equal  to  48  and  32.  This  degradation  of  the  array  gain  is 
consistent  with  the  estimates  of  the  reduction  of  the  processing  gain  obtained 
using  the  conventional  formulas  of  10  log  m  for  coherently  processing  m  sensors 
and  of  5  log  n  for  incoherently  averaging  over  n  estimates.  The  estimated 
processing  gain  reductions  are  2.7,  and  4.5  dB  for  processing  two  subapertures  with 
spatial  smoothing  lengths  equal  to  48  and  32  sensors.  In  the  case  of  the  high  input 
SNR,  the  array  gain  of  the  Minimum  Variance  processor  is  lower  as  a  result  of 
mismatch.  The  poor  performance  of  the  processor  may  also  be  due  to  a  signal 
cancellation  problem.  As  noted  earlier,  the  highest  level  arrival  actually  corresponds 
to  a  cluster  of  arrivals  very  close  together  which  is  hard  to  decorrelate  even  after 
heavy  spatial  smoothing.  The  reduction  of  the  array  gain  is  on  the  order  of  6  dB 
when  processing  two  subapertures. 


3.S.4.  Resolution  Study 

Using  the  MV  processor  results  in  some  processing  gain  degradation  which,  in 
some  cases,  is  compensated  for  by  a  better  resolution.  To  quantify  when  enhanced 
resolution  in  achieved  compared  to  the  conventional  processor,  a  simulation 
focussed  on  resolution  is  now  presented. 


3.5.4. 1.  Description  of  the  Study 

As  in  the  previous  sections,  the  Bai  .ett  processor  is  used  as  a  reference 
against  which  is  compared  the  MV  processor,  operating  on  the  full  aperture  (after 
spatial  smoothing  on  92,  78  and  64  sensors)  and  on  two  non-overlapping 
subapertures  (after  spatial  smoothing  on  48  and  32  sensors).  Resolution  for  each 
arrival  is  measured  here  through  a  relative  resolution  measure  with  respect  to  the 
conventional  beamformer.  The  relative  resolution  factor  for  a  particular  arrival  is 
defined  as 

BW(X)~  BW(  Bartlett  )  .  . 

BW(  Bartlett ) 

where  BW{X)  is  the  0.1  dB  beamwidth  of  the  considered  arrival  for  method  X.  When 
the  relative  resolution  factor  is  negative,  method  X  yields  better  resolution  than  the 
conventional  processor.  A  relative  resolution  factor  of  -1  means  that  the  resolution 
of  method  X  is  two  times  better  than  the  conventional  processor. 
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Since  it  was  found  that  there  is  no  advantage  of  using  curved  wavefront 
replica  vectors  for  the  900  m  aperture  array,  only  plane  wavefront  replicas  are  used 
in  this  resolution  study.  As  before,  4096  beams  between  -45°  and  45°  are  used, 
yielding  a  sampling  interval  in  angle  of  0.022°.  The  relative  beamwidth  is  computed 
from  the  0.1  dB  beamwidths  of  each  arrival  (interpolated  if  necessary)  and  plotted 
against  the  Bartlett  signal-plus-noise  to  noise  ratios  for  each  input  SNR  in  Figure 
3.19,  3.20  and  3.21. 


3. 5. 4. 2.  Results 


Full  Aperture  Processing 

The  MV  processor  produces  a  much  better  resolution  than  the  Bartlett 
processor  at  high  signal-to-noise  ratios.  This  is  illustrated  by  the  9  dB  input  SNR 
case  where  the  arrivals  with  output  signal-plus-noise  to  noise  ratios  larger  than  5 
dB  are  much  better  resolved  with  the  MV  processor,  whatever  the  spatial  smoothing 
length.  In  the  medium  input  SNR  case  (-11  dB),  the  resolution  has  decreased  to  the 
same  order  as  the  Bartlett  processor.  Using  a  spatial  smoothing  length  of  64  results 
in  a  degradation  of  resolution,  compared  to  full  aperture  conventional  processing. 
This  trend  is  confirmed  by  the  lower  input  SNR  case  (-18  dB),  where  a  larger  spatial 
smoothing  length  would  be  more  appropriate.  Therefore,  the  reduction  of  the 
efficient  aperture  length  because  of  spatial  smoothing  may  result  in  no  resolution 
improvement,  even  at  medium  input  SNRs. 


Two  Subaperture  Processing 

The  separation  of  the  array  into  two  subapertures  results  in  large  degradation 
of  resolution  as  shown  in  Figure  3.20.  The  MV  processor  operating  on  the 
subapertures  is  able  to  produce  a  better  resolution  than  the  conventional  processor 
only  for  arrivals  with  high  output  signal-plus-noise  to  noise  ratios.  At  low  input 
SNR,  the  MV  processor  produces  a  rather  poor  resolution. 


3.5.5.  Conclusions 

This  study  shows  that  full  aperture  MV  processing  can  provide  a  substantial 
improvement  in  resolution  compared  to  conventional  processing  at  high  (9  dB)  and 
medium  (-11  dB)  input  SNR,  but  that  it  does  not  for  the  low  (-18  dB)  input  SNR  case. 
The  improvement  in  resolution  for  large  level  arrivals  is  associated  with  a  large 
reduction  of  output  signal-plus-noise  to  noise  ratio,  as  a  result  of  mismatch  and 
signal  cancellation.  Slight  degradation  of  the  relative  peak  levels  for  low  level 
arrivals  was  observed  as  a  result  of  the  reduction  of  the  efficient  aperture  due  to 
spatial  smoothing. 

Subaperture  processing  on  two  subarrays  has  better  resolution  than 
conventional  beamforming  for  the  highest  level  arrivals  in  the  high  input  SNR  case. 
Splitting  the  array  in  two  halves  results  in  some  degradation  of  the  resolution  and 
also  some  loss  in  output  signal-plus-noise  to  noise  ratios,  with  a  half  dB  loss  for  the 
low  input  SNR  case,  a  few  dB  loss  for  the  medium  input  SNR  case,  and  up  to  a  6  dB 
loss  for  the  high  input  SNR  case. 
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Table  3.1 


Angle 

(deg) 

Loss 

(dB) 

0° 

0 

10° 

0.4 

20° 

3.20 

30° 

5.20 

40° 

6.80 

50° 

7.60 

60° 

8.20 

70° 

8.20 

80° 

8.40 

90° 

8.40 

:  Bottom  Loss  Table  (GSM  Province  Type  3  [Weinberg;  1985]). 
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Angle  (deg) 


Figure  3.1:  Snell’s  law  remapping. 
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Figure  3.5:  Conventional  beamformer  output.  The  curved  (solid  line)  and  plane  wave- 
front  (dotted  line)  angular  spectra  are  overlaid  with  the  GSM  eigenrays. 


i  ■  i  i 


-80  -60  -40  -20  0  20  40  60  80 

Angle  (deg) 


Figure  3.6:  Full  aperture  curved  wavefront  array  manifold. 
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Figure  3.7:  MV  processor  output  after  spatial  smoothing  on  64  sensors.  The  curved 
wavefront  MV  processor  output  (solid  line)  is  overlaid  with  the  curved  wavefront  con¬ 
ventional  processor  output  (dotted  line)  and  the  GSM  eigenrays. 


Figure  3.8:  Equivalent  Array  manifold  after  spatial  smoothing  on  64  sensors. 
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Figure  3.11:  Conventional  beamformer  output.  The  curved  (solid  line)  and  plane  wave- 
front  (dotted  line)  angular  spectra  are  overlaid  with  the  GSM  eigenrays. 
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Figure  3.12:  MV  processor  output  after  spatial  smoothing  on  64  sensors.  The  curved 
wavefront  MV  processor  output  (solid  line)  is  overlaid  with  the  curved  wavefront  con¬ 
ventional  processor  output  (dotted  line)  and  the  GSM  eigenrays. 
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subarrays 


Figure  3.15:  Peak  level  performance  study  for  a  -18  dB  input  SNR.  Left  panels:  full 
aperture  processing,  right  panels:  subaperture  processing.  Upper  panels:  eigensteering 
vectors,  lower  panels:  plane  wave  replicas. 


vectors,  lower  panels:  plane  wave  replicas. 


ire  3.17:  Peak  level  performance  study  for  a  9  dB  input  SNR.  Left  panels:  full  aper 
processing,  right  panels:  subaperture  processing.  Upper  panels:  eigensteering  vec 
lower  panels:  plane  wave  replicas. 
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Figure  3.18:  Array  Gain. 
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O  Smoothing  Length  -  64 
<5  Smoothing  Length  =  78 
A  Smoothing  Length  “  92 


Bartlett  S+N/N 


Bartlett  S+N/N 


Figure  3.19:  Resolution  performance  study  for  a  -18  dB  input  SNR.  Upper  panel:  full 
aperture  processing,  lower  panel:  subaperture  processing. 
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Figure  3.20:  Resolution  performance  study  for  a  -1 1  dB  input  SNR.  Upper  panel:  full 
aperture  processing,  lower  panel:  subaperture  processing. 
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Bartlett  S+N/N 


Figure  3.21:  Resolution  performance  study  for  a  9  dB  input  SNR.  Upper  panel:  full 
aperture  processing,  lower  panel:  subaperture  processing. 
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4.  Processing  Real  Data 


4.1.  Overview 

The  various  processing  techniques  proposed  in  Chapter  3  are  used  to  process 
real  data,  collected  at  sea  by  the  Marine  Physical  laboratory  digital  acoustic  array,  a 
900  m  long  aperture  with  120  sensors.  First,  the  data  set  is  presented  and  the 
acoustics  modelled  with  ray  and  normal  mode  models.  The  data  are  processed  by 
the  conventional  FFT  beamformer,  the  MV  processor  operating  on  the  full  aperture 
and  on  two  non-overlapping  subarrays  with  various  modified  spatial  smoothing 
lengths  and  plane  wavefront  replica  vectors.  Qualitative  results  are  presented  and 
analyzed  for  a  high  input  signal-to-noise  ratio  case  (11.8  dB).  A  quantitative  study 
completes  this  work  on  real  data  by  showing  that  in  both  high  and  low  input  SNRs 
(11.8  dB  and  -9.2  dB),  processing  gain  and  resolution  are  similar  to  the  patterns 
found  in  the  simulations  of  Chapter  3.  Full  aperture  MV  processing  provides  more 
resolution  but  less  processing  gain  than  conventional  processing.  Subaperture 
processing  does  not  perform  as  expected  because  of  the  great  multipath  variability 
of  the  particular  field  processed  here.  An  assumption  on  the  homogeneity  of  the 
pressure  field  has  been  implicitly  made  when  subaperture  processing  has  been 
formulated.  Experimental  results  using  these  techniques  indicate  that  real  pressure 
fields  can  be  highly  inhomogeneous. 


4.2.  Presentation  of  the  Data  Set 


4.2.1.  Introduction 

The  data  were  recorded  by  the  Marine  Physical  Laboratory  (MPL)  digital 
acoustic  array  during  the  September  1987  Vertical  Line  Array  (VLA)  Experiment  in 
the  North-East  Pacific  [Sotirin;  1989].  The  vertical  line  array,  900  m  long  and  filled 
with  120  hydrophones  with  7.5  m  interelement  spacing,  was  deployed  for  several 
weeks  from  the  Research  Platform  FLIP  in  a  three-point  moor  at  34°58.69'  N 
125°58.43'  W,  approximately  400  km  southwest  of  Monterey,  California.  The  VLA 
was  in  a  transponder  net  consisting  of  three  near-bottom  tranponders  to  enable 
acoustic  navigation  of  the  array  (Sotirin;  1989].  This  work  deals  with  the  200  Hz  CW 
tone  propagation  data  recorded  by  the  VLA  when  the  source  ship,  the  USNS 
NARRAGANSETT,  kept  a  fixed  station  89  nm  from  the  R/P  FLIP  on  27  September 
1987  or  Julian  Day  270  [Tran;  1989d].  The  top  sensor  of  the  VLA  was  at  a  nominal 
depth  of  400  m.  The  acoustic  source  was  deployed  at  20  m  depth  and  projected  a 
200  Hz  CW  tone  with  a  nominal  source  level  of  184  dB  ref  1  //Pa  at  1  m  for  18 
minutes  from  21:57  GMT.  The  data  were  recorded  on  Tape  918  (which  starts  at 
21:58:27  GMT). 


4.2.2.  Description  of  the  VLA  System 

The  characteristics  of  the  MPL  digital  array  are  summarized  here,  more  detailed 
information  can  be  found  in  [Sotirin;  1988].  The  vertical  line  array  has  12  sections, 
each  section  contains  10  acoustic  transducers  and  a  microprocessor-based 
processor  which  samples  and  telemeters  the  data  to  topside  on  a  serial  link.  The 
acoustic  channels  are  numbered  from  bottom  to  top  and  a  processor  samples  the  5 
transducers  immediately  below  and  above.  The  signal  in  a  given  channel  is  the 
output  of  a  set  of  two  Benthos  Aquadyne  AQ-1  hydrophones  wired  in  series.  The 
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resulting  acoustic  transducer  has  a  sensitivity  of  -197  dB////Pa.  The  transducer 
output,  then,  is  amplified  with  a  pre-amplification  stage  with  40  dB  gain  and  filtered 
between  10  and  220  Hz  (with  a  nominal  bandpass  gain  of  0  dB).  A  differential  line 
driver  transmits  the  signal  to  the  processor  producing  a  6  dB  gain.  In  the 
processor,  the  analog  signal  is  amplified  by  a  software  controlled  variable  gain 
(selectable  from  0  to  60  dB).  The  10  acoustic  channels  are  multiplexed  and  sampled 
by  a  12  bit  bipolar  A/D  converter  with  a  ±10  V  reference  voltage.  The  sampling  rate 
is  500  Hz.  The  two  channels  immediately  adjacent  to  a  processor  are  low-pass 
filtered  in  the  processor  to  allow  the  12  kHz  acoustic  navigation  signals  to  enter 
into  the  detection  circuitry  in  the  telemery  module. 


4.2.3.  Environmental  Data 

Environmental  data  were  collected  during  the  experiment  on  the  R/P  FLIP. 

Swell  heights  on  Julian  Day  270  were  visually  measured  between  5  and  10  feet.  A 
northern  wind  was  blowing  with  speeds  between  20  and  25  knts.  A  CTD  was 
collected  in  the  vicinity  of  the  R/P  FLIP  on  Julian  Day  267,  up  to  a  depth  of  3885  m. 
Another  CTD,  collected  at  37°04.19'  N  134°46.75'  W  on  Julian  Day  264,  was  used  to 
extend  the  CTD  at  the  R/P  FLIP  down  to  the  nominal  bottom  depth,  4667  m.  The 
sound  speed  profile,  based  on  the  CTD,  was  derived  using  the  UNESCO  equations 
relating  conductivity  to  practical  salinity  and  sound  speed  to  temperature,  salinity 
and  pressure  [Fofonoff;  1983].  XBT  casts  were  collected  at  the  R/P  FLIP  before  and 
after  the  fixed  station  transmissions  on  Julian  Day  270  at  13:00  GMT  (a  400  m  XBT) 
and  on  Julian  Day  271  at  03:00  GMT  (a  750  m  XBT).  These  XBT  data  were  used  in  the 
state  equation  of  [Mackenzie:  1981]  to  produce  sound  speed  profiles.  This 
information  allows  the  synthesis  of  a  composite  sound  speed  profile  that  will  be 
used  in  the  acoustic  modeling.  A  desampled  version  of  the  400  m  XBT  sound  speed 
profile  is  used  from  0  to  250  m.  Then,  a  desampled  version  of  the  750  m  XBT  is 
used  down  to  700  m  where  it  blends  with  the  CTD  sound  speed  profile.  A 
desampled  version  of  the  CTD  sound  speed  profile  is  then  used  down  to  the 
bottom.  The  model  sound  speed  profile  is  plotted  in  Figure  4.1  for  the  whole  water 
column  and  for  the  top  750  m.  The  two  XBT  sound  speed  profiles  are  also  plotted 
in  dashed  lines.  A  30  m  deep  mixed  layer  was  observed  at  the  R/P  FLIP,  followed  by 
a  shallow  and  weak  duct  above  100  m,  which  evolves  with  time,  as  shown  by  the 
two  XBT  sound  speed  profiles  taken  only  14  hours  apart.  The  environment  has  a 
deep  sound  axis  around  550  m. 


4.2.4.  Data  Quality  Checks 

Preliminary  data  quality  checks  included  time  series  plots,  power  spectra  for 
selected  array  channels  and  conventional  FFT  beamformer  outputs.  The  power 
spectra,  computed  between  DC  and  250  Hz  and  plotted  in  Figure  4.2,  are  obtained 
by  processing  2  minutes  and  11  seconds  of  data  They  are  derived  from  the 
incoherent  addition  of  15,  50  %  overlapped,  8192-point  FFTs  (61  mHz  bin  width).  A 
Kaiser-Bessel  window  with  a  parameter  of  2.5,  yielding  a  sidelobe  level  of  -  57  dB 
[Harris:  1978],  weights  the  data  prior  to  each  FFT.  Power  values  are  calibrated  in  dB 
re  1  ^Pa,  to  yield  the  power  for  sinusoidal  signals.  The  90  %  confidence  interval  in 
these  spectra  is  +2/-1.6  dB.  The  power  spectra  illustrate  the  good  quality  of  the 
recorded  data.  The  200  Hz  signal  appears  as  a  strong  line  component.  The  60  Hz 
line  and  its  harmonics  (120  Hz,  180  Hz  and  240  Hz)  are  the  result  of  power  line 
contamination.  Some  variability  in  the  noise  levels  at  high  frequency  is  observed, 
the  channels  have  different  self  noise  characteristics. 
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The  distributions  of  power  in  the  vertical,  plotted  in  Figure  4.3,  correspond  to 
the  200  Hz  signal  bins  7373  and  7374  (199.951  Hz  and  200.134  Hz),  and  the  noise 
bins  7368  and  7380  (199.646  Hz  and  200.378  Hz).  They  are  computed  by 
incoherently  averaging  thirty  50  %  overlapped  8192-point  FFTs  obtained  by 
processing  two  data  segments  of  2  minutes  and  11  seconds.  As  before,  the  data  are 
weighted  by  a  Kaiser-Bessel  window  with  a  parameter  of  2.5  prior  to  the  FFT.  The 
90  %  confidence  interval  in  Figure  4.3  is  +  1.4/-1.2dB.  The  maximum  signal  power 
varies  between  80  and  90  dB  re  1  /iPa  which  is  consistent  with  a  source  level  of  184 
dB///iPa  at  1  m  that  suffers  spherical  transmission  loss  over  a  range  on  the  order  of 
a  water  depth  (5000  m)  and  then  cylindrical  loss  over  the  remaining  range.  The 
distribution  of  noise  power  across  the  array  shows  a  somewhat  large  variability, 
with  peaks  5  to  15  dB  above  the  noise  floor  between  60  and  65  dB  (re  1  pi  Pa).  The 
normalization  factor  to  pass  from  the  sinusoidal  component  normalization  to  a 
noise  spectral  density  normalization  is  on  the  order  of  10  dB  so  the  noise  floor 
level  at  200  Hz  is  between  70  and  75  dB//p<Pa /^Hz.  Such  levels  are  consistent  with 
the  common  noise  levels  summarized  in  [Urick;  1986]  for  the  wind  speeds  observed 
at  the  R/P  FLIP  (20  to  25  knts).  The  peaks  in  the  distribution  of  noise  power  in  the 
vertical  are  not  the  result  of  a  statistical  artifact  given  the  large  number  of  averages 
in  the  estimator.  As  mentioned  earlier,  they  indicate  channels  with  high  self  noise 
characteristics.  The  peaks  generally  correspond  to  the  channels  which  are  physically 
the  closest  or  the  farthest  to  the  processors.  The  sensors  which  are  closest  to  the 
processors  are  used  for  the  array  navigation  and  are  known  to  be  more  noisy.  The 
sensors  farthest  from  the  processor  in  each  75  m  section  suffer  the  most  cross-talk 
with  other  telemetry  and  transmission  lines. 

The  FFT  beamformer  processes  a  signal  and  a  noise  FFT  bins  (the  signal  bin  is 
7373  and  the  noise  bin  is  7368)  to  obtain  preliminary  information  on  the  vertical 
arrival  structure.  The  array  is  cut  for  100  Hz  with  an  interelement  spacing  of  7.5  m 
and  is  therefore  spatially  aliased  at  the  operating  frequency  of  200  Hz.  The 
beamforming  results  are  thus  plotted  between  ±29.5°  with  respect  to  the  horizontal 
in  Figure  4.4.  Since  the  array  is  numbered  from  bottom  to  top,  positive  angles 
correspond  to  downgoing  sound  or  uplooking  beams.  The  noise  arrival  structure  is 
uniform  with  angle.  The  90  %  confidence  interval  is  2/- 1.6  dB  except  for  the 
spectra  referenced  at  22:00:38,  22:13:48  and  22:19:12  GMT  where  it  is  2.5/-1.8,  3.3/- 
2.3  and  1.8/- 1.44  dB,  respectively.  The  arrival  structure  is  stationary  since  the  ship 
was  keeping  fixed  station.  Only  steep  arrivals  around  ±15°  are  observed,  they 
correspond  to  refracted  surface-reflected  (RSR)  paths. 


4.2.5.  Acoustic  Modeling 


Modeling  with  the  GSM 

The  CONGRATS  ray  tracing  program  of  the  Generic  Sonar  Model  [Weinberg; 
1985]  is  utilized  to  model  the  sound  field.  The  model  sound  speed  profile  of 
Section  4.2.3  is  used  and  the  bottom  corresponds  to  the  GSM  Province  Type  3  (c.f. 
Table  3.1).  The  ray  traces  in  Figure  4.5  for  a  20  m  deep  source  shows  that 
convergence  zone  propagation  takes  place  with  refracted  surface-reflected  rays.  At 
165  km  from  the  source,  the  array  receives  sound  that  travelled  over  3  convergence 
zones.  Although  the  nominal  source  level  is  184  dB  re  1  piPa  at  1  m,  the  calibrated 
received  levels  obtained  from  the  experimental  data  indicate  that  the  source  level  is 
actually  closer  to  174  dB  re  1  piPa  at  1  m  [Tran;  1989d].  This  latter  source  level  will 
be  used  in  the  following.  Using  a  trial  and  error  procedure,  a  best  fit  between  the 
experimental  data  and  the  modeling  is  obtained  for  a  source  range  of  162  km.  The 
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distribution  of  power  computed  by  GSM.  plotted  in  Figure  4.6,  exhibits  a  strong 
interference  pattern  similar  to  the  data.  The  GSM  eigenray  angles  at  the  sound  axis 
(550  m)  are  plotted  in  Figure  4.7  at  162  km.  Once  again,  a  good  correspondance  with 
the  FFT  beamformer  output  is  observed. 


Modeling  with  ATLAS 

The  modeling  is  now  performed  with  the  ATLAS  normal  mode  model  [Gordon; 
1984].  The  environmental  parameters  are  the  same  as  with  the  GSM.  The  power  as  a 
function  of  range,  computed  by  the  ATLAS  model,  is  plotted  in  Figure  4.8  for  a 
receiver  a,.  the  same  depth  as  the  200  Hz  source  (i.e.  20  m).  This  shallow  depth 
enhances  the  appearance  of  the  convergence  zones.  The  sound  propagates  over 
three  convergence  zones  from  the  source  to  the  receiving  array,  assumed  at  165  km 
in  the  trailing  edge  of  the  third  convergence  zone.  The  model  distribution  of  power 
across  the  array  is  plotted  in  Figure  4.9  for  a  source  at  161.4  km  from  the  VIA,  and 
the  modeling  is  excellent  as  with  the  GSM.  The  ATLAS  normal  mode  program 
produces  the  complex  wavefield  across  the  vertical  array  which  can  be  beamformed 
to  yield  the  model  vertical  arrival  structure.  The  synthetic  wavefield  for  a  165  km 
range  is  beamformed  as  with  the  experimental  data  and  plotted  in  Figure  4.10.  A 
very  good  match  of  the  arrival  structures  is  observed.  The  pair  of  arrivals  above 
10°  corresponds  to  downgoing  sound,  which  is  consistent  with  having  the  receiving 
array  in  the  trailing  edge  of  the  third  convergence  zone. 


4.3.  Processing  Real  Data 

The  three  processing  approaches  proposed  in  Chapter  3  are  now  used  on  these 
data  collected  at  sea.  Only  plane  wavefront  replica  vectors  will  be  used  with  the 
Minimum  Variance  processor  operating  on  the  full  aperture  and  on  two  non¬ 
overlapping  subapertures.  4096  beams  are  used  between  -29.5°  and  29.5°.  The  first 
step  of  the  processing  is  to  do  a  time  to  frequency  domain  transformation.  The  time 
series  of  each  channel  are  divided  into  50  %  overlapped  data  blocks  called  time 
snapshots  and  Fourier  transformed  using  a  Fast  Fourier  Transform  (FFT).  The  FFT  is 
performed  here  with  a  Kaiser-Bessel  window  with  a  parameter  of  2.5  (-57  dB 
sidelobe  level).  The  FFT  length  determines  the  input  signal-to-noise  ratio  (SNR),  i.e. 
the  signal-to-noise  ratio  at  the  input  of  the  beamformer.  The  preliminary  analysis 
presented  in  Section  4.2  showed  that  the  input  SNR  is  on  the  order  of  17  dB  for  a 
8192-point  FFT.  Since  coherent  processing  with  a  FFT  length  N fft  yields  a  processing 
gain  of  10  log  N fft,  the  FFT  length  is  reduced  to  2048  and  16  points  so  that  the 
input  SNRs  are  on  the  order  of  11  and  -10  dB,  respectively,  and  similar  to  those  of 
Chapter  3.  For  a  2048-point  FFT,  the  bin  width  is  244  mHz  and  the  200  Hz  signal 
lays  in  bin  1844  (natural  frequency  of  200.2  Hz).  For  a  16-point  FFT,  the  bin  width  is 
31.25  Hz  and  the  200  Hz  signal  is  in  bin  15  (center  frequency  of  218.75  Hz). 

The  conventional  beamformer  does  an  across  FFT  on  the  selected  frequency 
bin  with  a  Kaiser-Bessel  window  with  a  parameter  of  1.5  (to  yield  a  35  dB  sidelobe 
level  rejection).  The  dyad  vector  formed  by  the  signal  frequency  bin  across  the 
array  is  Fast  Fourier  Transformed  to  produce  a  periodogram.  The  periodograms  are 
averaged  for  a  number  of  snapshots  to  obtain  a  statistically  stable  angular 
spectrum.  Here,  fifteen  averages  are  performed  and  the  90  %  confidence  interval  is 
+2/- 1.6  dB. 

The  Minimum  Variance  processor  requires  the  estimation  of  the  array 
covariance  matrix.  The  covariance  matrix  is  computed  by  averaging  over  time 
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snapshots  the  dyad  vectors  computed  by  the  time  to  frequency  transform.  At  least 
as  many  dyads  as  sensors  must  be  averaged  to  ensure  full  rank.  The  MV  processor, 
when  it  processes  the  full  aperture,  requires  at  least  120  averages  in  the  estimate 
the  array  covariance  matrix.  Here,  126  snapshots  are  used.  When  the  MV  processor 
beamforms  the  two  60-sensor  subapertures,  at  least  60  averages  are  required  to 
estimate  the  array  covariance  matrices  and  63  snapshots  are  used  here.  In  all  cases 
the  mairix  is  readily  invertible  after  stabilization  with  a  fraction  of  noise  equal  to 
10-6.  The  time  duration  of  the  data  blocks  processed  by  the  three  methods  is 
summarized  in  Table  4.1. 

Since  the  array  receives  multipath  arrivals,  spatial  smoothing  is  performed 
prior  to  using  the  MV  processor.  When  the  full  aperture  is  processed,  three  spatial 
smoothing  lengths  are  used:  92,  78  and  64.  When  the  two  60-sensor  subapertures 
are  processed,  two  spatial  smoothing  lengths  are  used:  48  and  32. 

The  actual  count  of  operations  for  each  processing  step  of  the  Minimum 
Variance  processor  is  examined  in  [Tran:  1989d].  Given  the  large  number  of  beams, 
the  most  computationally  intensive  operation  is  the  quadratic  form.  Neglecting  the 
additional  load  required  to  map  one  of  the  subaperture  angular  spectra  using 
Snell’s  law  before  incoherently  averaging,  subaperture  processing  after  spatial 
smoothing  on  48  sensors  was  found  to  reduce  the  computational  load  by  25  % 
compared  to  full  aperture  MV  processing  after  spatial  smoothing  on  78  sensors. 


Method 

Number 
of  FFTS 

Data  block  I 
2048-point  FFT 

ength  for 
16-point  FFT 

Conventional  FFT 
Beamformer 

15 

32.8  sec. 

0.256  sec. 

MV  processor 
on  full  aperture 

126 

4  min.  20  sec. 

2  sec. 

MV  processor 
on  two  subapertures 

63 

2  min.  1 1  sec. 

1  sec. 

Table  4.1:  Length  of  the  data  blocks 

4.4.  Qualitative  Study  of  the  MV  Processing  Results:  1 1  dB  Input  SNR 


4.4.1.  Conventional  Processing 

The  conventional  FFT  beamformer  output  corresponding  to  the  FFT  bin  1844 
(center  frequency  of  200.2  Hz),  is  plotted  in  Figure  4.1 1.  The  plane  wave  replica 
vectors  are  referenced  to  the  sound  axis  sound  speed  (1479.1  m/s).  There  are  two 
main  downgoing  arrivals  in  the  spectrum.  The  relative  signal  level  in  Figure  4.11  is 
only  on  the  order  of  25  dB  and  sidelobe  leakage  is  experienced,  although  a  Kaiser- 
Bessel  window  with  a  parameter  of  1.5  with  a  35  dB  sidelobe  rejection  is  used.  The 
array  beampattern  is  degraded  because  of  phase  and  gain  errors  from  channel  to 
channel  [Sotirin;  1989]. 


4.4.2.  Full  Aperture  MV  Processing 

Full  aperture  MV  processing  is  performed  on  bin  number  1844  after  spatial 
smoothing  on  the  array  covariance  matrix.  Plane  wavefront  replica  vectors  are  used 
and  referenced  to  the  sound  axis  sound  speed.  When  the  spatial  smoothing  length 
is  64,  a  large  number  of  averages  is  performed  to  decorrelate  the  multipath  arrivals. 
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One  observes  that  the  MV  estimate  in  Figure  4.12  still  suffers  large  signal 
cancellation.  The  arrivals  between  10°  and  15°  are  close  together  and  difficult  to 
decorrelate.  Since  plane  wavefront  replica  vectors  are  used,  modified  spatial 
smoothing  should  produce  better  results  and  be  used.  The  MV  processor  output, 
obtained  after  modified  spatial  smoothing  is  plotted  in  Figure  4.13.  The  power 
levels  are  now  closer  (at  most  a  few  dB  lower)  to  the  ones  produced  by  the  FFT 
beamformer.  Therefore,  modified  spatial  smoothing  will  be  used  in  the  following. 
When  the  modified  spatial  smoothing  lengths  are  respectively  78  and  92,  resolution 
is  enhanced  at  the  cost  of  more  signal  cancellation  as  shown  in  Figure  4.14  and  4.15. 
The  arrival  around  12°  as  well  as  the  one  around  14°  actually  corresponds  to  a 
cluster  of  at  least  tv/o  arrivals  very  close  in  angle.  A  subjective  inspection  of  the  full 
aperture  MV  processing  results  suggests  that  a  modified  spatial  smoothing  length 
of  78  is  optimal. 


4.4.3.  Subaperture  MV  Processing 


4.4.3. 1.  Experimental  Results 

Subaperture  MV  processing  is  performed  on  two  non-overlapping  60-sensor 
subarrays.  The  upper  half  of  the  array  extends  from  400  to  842.5  m,  and  is 
referenced  to  a  harmonic  sound  speed  on  the  order  of  1479.1  m/s  (the  sound  axis 
sound  speed).  The  lower  half  of  the  array  extends  from  850  m  to  1292.5  m  and  is 
referenced  to  a  harmonic  sound  speed  of  1481.4  m/s.  The  two  subarrays  are 
processed  separately  using  their  respective  plane  wavefront  replica  vectors.  A 
modified  spatial  smoothing  length  of  48  is  used.  Since  the  lower  aperture  is 
referenced  to  a  lower  sound  speed,  its  angular  spectrum  is  remapped  from  a 
reference  sound  speed  of  1481.4  to  1479.1  m/s  using  Snell’s  law.  The  upper 
subarray  angular  spectrum  and  the  remapped  lower  subarray  angular  spectrum  are 
plotted  in  Figure  4.16  and  4.17,  respectively.  The  composite  angular  spectrum  is 
obtained  by  incoherently  averaging  the  two  subaperture  spectra  and  is  plotted  in 
Figure  4.18.  Three  arrivals  are  resolved  between  10°  and  15°  instead  of  the  two 
groups  of  closely  separated  arrivals  detected  by  the  full  aperture  MV  processor. 
These  results  are  unexpected,  the  upper  and  lower  subaperture  angular  spectra  are 
not  aligned,  even  using  Snell’s  law.  Subaperture  and  full  aperture  processing 
produce  different  arrival  structures.  The  upper  half  of  the  array  receives  a  high 
power  level  arrival  around  14°  and  a  lower  level  one  around  12.7°.  The  lower  half 
receives  two  arrivals  around  12.3°  and  12.7°. 

Using  a  modified  spatial  smoothing  length  of  32  produces  similar  results, 
although  the  estimate  of  the  angular  spectrum  is  of  poor  quality  as  a  result  of  the 
short  smoothing  length.  Between  10°  and  15°,  the  upper  subaperture  detects 
steeper  arrivals  than  the  lower  aperture  as  with  smoothing  on  48  sensors.  The 
resulting  composite  spectrum  exhibits  a  peak  around  -15°  and  a  broad 
conglomerate  of  peaks  between  10°  and  15°.  A  subjective  look  at  these  results 
suggests  that  using  a  spatial  smoothing  length  of  32  sensors  is  not  adequate.  The 
emphasis  is  now  put  on  the  acoustic  modeling,  in  an  attempt  to  understand  the 
behaviour  of  the  subaperture  processor. 


4. 4. 3. 2.  Acoustic  Modeling 


Using  GSM 

The  preliminary  acoustic  modeling,  presented  in  Section  4.2.5,  is 
complemented  by  a  more  detailed  study  of  the  eigenrays  between  10°  and  15°  for 
receiver  depths  between  200  to  1600  m.  As  before,  the  source  is  at  20  m  depth  and 
a  162  km  range.  All  possible  eigenrays  within  the  [10M5°]  interval  are  plotted  as  a 
function  of  depth  in  Figure  4.19.  Because  of  the  variations  of  sound  speed  with 
depth,  the  eigenrays  at  each  receiver  depth  arrive  at  different  angles  which  stay 
within  the  studied  angular  interval.  Figure  4.19  shows  that  four  dominant  arrivals 
impinge  on  the  array,  they  correspond  to  two  pairs  of  closely  separated  arrivals 
which  are  not  resolved.  At  400  m  depth,  the  first  arrival  pair  has  arrival  angles 
equal  to  1 1.21°  and  1 1.24°  and  power  levels  equal  to  67.  dB  and  70.  dB.  The  two 
other  arrivals  are  steeper  at  13.3°  and  13.4°  and  have  levels  of  75.  and  76.  dB.  The 
steeper  pair  in  the  [10°,15°]  interval  dominates  the  upper  half  of  the  array.  Over 
the  lower  half,  both  pairs  are  present. 

Ray  traces  for  take-off  angles  between  -15°  and  15°  with  a  1°  increment, 
plotted  in  range  between  160  and  164  km  and  in  depth  between  200  to  1600  m  in 
Figure  4.20.  show  the  same  phenomena.  An  aperture  deployed  at  162  km  receives 
steeper  rays  between  400  to  842.5  m  than  in  the  850  to  1292.5  m  depth  interval. 
This  modeling  with  the  GSM  indicates  that  the  variability  observed  in  the 
experimental  arrival  structure  is  real. 


Modeling  with  ATLAS 

The  ATLAS  normal  mode  model  is  now  used  to  check  these  preliminary 
conclusions.  The  environmental  model  is  the  same  as  in  the  preliminary  acoustic 
modeling.  First,  the  implementation  of  the  Snell’s  law  remapping  operation  is 
investigated.  A  canonical  ATLAS  pressure  field  corresponding  to  a  single  mode  is 
beamformed  with  the  two  subaperture  MV  processor.  The  normalized  modal 
eigenfunctions,  gray-level  displayed  in  Figure  4.21  for  the  whole  water  column  and 
for  mode  number  3  to  303,  summarizes  a  great  deal  of  the  acoustics.  It  shows  which 
modes  are  excited  by  the  20  m  source  (modes  around  mode  number  147)  and  which 
ones  are  bottom  interacting.  The  field  corresponding  to  mode  number  147  is  used 
to  create  a  simulated  array  covariance  matrix 

R  =  a}  d  dH  +  a2n  I  (4.1) 

where  the  d  is  the  normalized  ATLAS  pressure  field  such  that  dH  d  =  1,  and  I  is  the 
identity  matrix.  An  input  signal-to-noise  ratio  of  10  dB  and  a  noise  level  of  50  dB 
are  assumed  so  that  the  simulated  noise  power  is  <xf  =  120.  10s  and  the  simulated 
signal  level  is  erf  =  14.4  109.  Since  a  single  mode  is  excited,  the  pressure  field  has  a 
constant  phase  modulo  ir.  Each  mode  corresponds  to  an  upgoing  ray  and  a 
downgoing  ray  [Boyles;  1984]  as  observed  in  Figure  4.22  where  are  plotted  the 
upper  and  lower  subaperture  spectra  produced  by  the  MV  processor  after  spatial 
smoothing  on  a  48  sensors.  Since  the  lower  aperture  is  referenced  to  a  lower  sound 
speed,  its  arrival  angles  are  shallower.  The  Snell’s  law  remapping  aligns  perfectly 
the  peaks  of  the  upper  and  lower  subaperture  spectra,  as  shown  in  Figure  4.23 
where  the  two  subaperture  spectra  are  overlaid. 

The  transmission  loss  from  10  to  1600  m  in  depth  and  from  163  to  167  km  in 
range  is  gray-level  displayed  between  40  and  65  dB  in  Figure  4.24.  Dark  areas 
correspond  to  regions  of  strong  attenuation  and  light  colored  areas  correspond  to 
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region  of  weak  attenuation.  The  best  fit  between  the  ATLAS  angular  spectrum  and 
the  data  was  found  for  a  range  of  165  km.  There,  the  array  is  in  the  trailing  edge  of 
the  third  convergence  zone.  White  lineations,  which  are  downgoing  with  range, 
correspond  to  the  downgoing  ray  paths  of  the  GSM.  The  400-850  m  depth  interval 
is  insonified  by  a  sharp  and  steep  light  colored  line  while  the  850-1300  m  interval 
is  insonified  by  a  broad  light-colored  area  which  has  smaller  slope  with  range.  The 
ATLAS  field  plot  is  in  very  good  agreement  with  the  GSM  ray  traces:  steeper 
arrivals  dominate  the  upper  half  of  the  array.  To  confirm  further  the  agreement 
between  the  modeling  and  the  data,  the  field  produced  by  the  ATLAS  model  across 
the  array  at  165  km  is  beamformed  by  the  subaperture  MV  processor  (with  a 
modified  spatial  smoothing  length  of  48).  The  simulated  field  covariance  matrix  is 
given  by  Equation  (4.1),  using  the  same  noise  and  signal  power.  The  simulated 
upper  and  lower  subapertures  are  processed  separately  and  their  angular  spectra 
plotted  in  Figure  4.25.  The  composite  spectrum  is  plotted  in  Figure  4.26.  In  the 
[10°,15°]  interval,  the  phenomenon  observed  in  the  data  is  also  present  in  the 
simulated  data.  Steeper  arrivals  dominate  in  the  upper  half  of  the  array  while  both 
shallower  and  steeper  arrivals  are  present  in  the  lower  half  of  the  array. 


4.4.4.  Conclusions 

Full  aperture  processing  produces  a  large  increase  of  resolution,  as  expected  in 
this  high  signal-to-noise  ratio.  The  cost  is  a  reduction  of  the  processing  gain.  The 
arrivals  between  10°  and  15°  are  composed  of  a  cluster  of  arrivals  very  close 
together  and  difficult  to  decorrelate  with  spatial  smoothing  and  modified  spatial 
smoothing  is  the  prefered  smoothing  procedure.  Subaperture  processing  does  not 
behave  as  expected.  An  assumption  on  the  homogeneity  of  the  field  has  been 
implicitely  made  when  subaperture  processing  was  proposed.  A  great  variability  in 
the  multipath  structure  was  found  in  the  experimental  data  and  was  succesfully 
modeled  using  the  GSM  and  the  ATLAS  models.  The  angular  spectra,  computed 
across  the  upper  and  lower  halves  of  the  array,  can  be  so  different  that  averaging 
does  not  enhance  the  peaks.  The  composite  subaperture  angular  spectrum,  for  a 
spatial  smoothing  length  of  48,  have  three  peaks  instead  of  two  between  10°  and 
15°.  Two  distinct  arrivals  are  reported  around  12°  where  an  unresolved  cluster  of 
arrivals  is  reported  by  full  aperture  processing.  On  the  other  hand,  subaperture 
processing  resolves  a  single  arrival  around  14°  while  full  aperture  processing  MV 
processing  indicates  very  closely  separated  arrivals. 


4.5.  Quantitative  Study 

Performance  aspects  are  now  studied  as  in  the  simulations  presented  in 
Chapter  3,  for  the  11  and  -9  dB  input  signal-to-noise  ratios.  Two  performance 
indices  are  examined:  for  each  arrival,  the  signal-plus-noise  to  noise  ratio  and  a 
relative  resolution  factor  which  is  the  relative  difference  in  0.1  dB  beamwidth 
between  the  MV  and  the  Bartlett  angular  spectra.  Here,  the  goal  of  the  study  is  to 
relate  the  performance  measures  on  real  data  to  the  performance  patterns  outlined 
through  simulations  in  Chapter  3.  Since  the  signal  frequencies,  arrival  structures, 
and  the  array  geometries  are  different  in  the  data  and  in  the  simulations,  only  a 
partial  agreement  is  expected. 
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4.5.1.  11  dB  Input  SNR 

The  three  arrivals  around  -15°,  12°  and  14°,  modelled  with  GSM  and  ATLAS,  are 
considered  here.  Their  peak  power  is  easily  measured  from  the  angular  spectra 
presented  in  the  previous  Sections.  In  the  case  of  subaperture  MV  processing  with  a 
modified  spatial  smoothing  length  of  48,  three  arrivals  are  detected  instead  of  two 
between  10°  and  15°.  Based  on  the  acoustic  modeling  described  in  Section  4.43.2, 
the  pair  of  peaks  around  12°  is  considered  to  correspond  to  the  12°  cluster  of 
arrivals,  and  its  higher  level  peak  is  selected  as  the  12°  arrival  and  used  to  measure 
signal-plus-noise  level  and  0.1  dB  beamwidth. 

The  estimation  of  the  output  noise  power  requires  a  200  Hz  signal  free  noise 
data.  Signal  which  is,  in  frequency,  five  bins  away  from  the  200  Hz  signal  bin 
(number  1844)  is  used  to  compute  an  estimate  of  the  noise  covariance  matrix.  The 
plane  wavefront  replica  vectors  are  computed  to  match  the  frequency  of  the  noise 
bin  (number  1839)  or  198.73  Hz.  The  beamformed  noise  field  is  evenly  distributed 
in  angle  (spatially  white)  which  is  consistent  with  the  mid-frequency  considered 
here,  a  spatially  aliased  array  and  measured  wind  speed  around  25  knts.  The 
signal-plus-noise  to  noise  ratios  for  the  three  detected  arrivals  are  compared  to  the 
corresponding  signal-plus-noise  to  noise  ratios  of  the  FFT  beamformer.  Their 
difference  is  plotted  in  Figure  4.27  together  with  the  simulations  results  of  [Tran; 
1989b],  where  the  input  SNR  was  10  dB.  The  simulation  data  are  plotted  with  the 
small  size  symbols  which  corresponds  to  a  given  spatial  smoothing  length.  Original 
spatial  smoothing  was  performed  in  the  simulations  while  modified  spatial 
smoothing  is  used  with  the  experimental  data.  Although  signal  frequencies,  arrival 
structure  and  array  geometries  are  different,  good  consistency  is  observed  between 
the  experimental  and  the  simulation  results.  The  relative  peak  levels  have  the  same 
order  as  in  the  simulations.  For  full  aperture  MV  processing,  the  drop  in  signal- 
plus-noise  to  noise  ratio,  compared  to  conventional  beamforming,  is  at  most  6  dB. 
The  experimental  results  suggest  that  there  is  still  some  signal  cancellation  for  the 
high  level  arrivals  since  the  relative  peak  levels  increase  when  the  spatial 
smoothing  length  increases.  The  subaperture  MV  processing  results  are  in  good 
agreement  with  the  simulation,  with  slightly  worse  performance  for  the  low  power 
-15°  arrival.  The  relative  peak  levels  of  the  highest  power  arrivals  are  consistent 
with  the  simulation  results,  differences  may  be  due  to  the  multipath  variability 
observed  in  the  data. 

The  relative  resolution  factor  is  derived  from  the  0.1  dB  beamwidth  measured 
in  the  angular  spectra  using  Equation  (3.6).  It  is  plotted  against  the  corresponding 
signal-plus-noise  to  noise  ratios  at  the  output  the  conventional  FFT  beamformer  in 
Figure  4.27.  The  simulation  results  (corresponding  to  the  10  dB  input  SNR)  are 
overlaid  with  smaller  symbols.  Once  again,  one  observes  a  good  consistency 
between  the  experimental  results  and  the  data,  with  a  slighty  worse  resolution  in 
the  experimental  results.  Still,  full  aperture  MV  processing  produces  an 
improvement  in  resolution.  Subaperture  processing  produces  mixed  results  as  a 
result  of  the  multipath  variability. 


4.5.2.  -10  dB  Input  SNR 

Only  two  arrivals  are  detected  by  the  various  beamformers  between  10°  and 
15°  in  this  lower  input  SNR  case.  Since  the  FFT  used  here  is  only  16  point  long,  the 
FFT  binwidth  is  30  Hz  wide  and  the  noise  power  cannot  be  calculated  as  in  the 
previous  Section.  The  200  Hz  transmissions  lasted  for  18  minutes  and  stopped  at 
22:17:01  GMT.  The  data  at  22:20:17  GMT,  free  from  the  200  Hz  CW  tone,  are  used  to 
estimate  the  noise  angular  spectra.  The  data  are  Fourier  transformed  with  a  16 
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point  FFT  and  bin  number  15  is  beamformed.  as  with  the  200  Hz  signal.  Although 
noise  data  were  recorded  18  minutes  after  the  signal  data,  this  estimate  of  the  noise 
power  should  have  limited  bias  and  be  reliable. 

The  signal-plus-noise  to  noise  ratios  are  computed  and  plotted  in  Figure  4.28 
against  the  conventional  FFT  beamformer  signal-plus-noise  to  noise  ratios.  The 
experimental  results  are  close  to  the  simulation  results  in  the  case  of  full  aperture 
MV  processing,  with  very  limited  processing  gain  degradation.  The  two  subaperture 
MV  processing  results  produces  slightly  better  performance  than  conventional 
beamforming  on  the  full  aperture  for  the  arrival  with  relative  peak  level  equal  to  4.5 
d3.  When  a  subsegment  length  of  32  is  used,  only  the  12°  arrival  is  resolved. 

The  relative  resolution  factor  is  computed  as  before  and  plotted  in  Figure  4.28. 
Once  again,  the  experimental  results  are  in  good  agreement  with  the  simulation 
results  (plotted  in  small  size).  Full  aperture  processing  produces  an  increase  in 
resolution  for  the  6  dB  signal-plus-noise  to  noise  ratio  arrival.  The  two  subaperture 
MV  processor  operating  on  the  data  has  resolution  characteristics  similar  to  the 
simulation  results.  But  the  highest  arrival  is  slightly  better  resolved  than  with  the 
conventional  FFT  beamformer. 


4.6.  Array  Gain  Calculation 

The  array  gain  obtained  by  the  various  beamformers  can  be  computed 
assuming  that  the  detection  is  based  on  the  highest  level  arrival.  The  input  and 
output  signal-to-noise  ratios  (SNR)  are  calculated  from  the  signal-plus-noise  and 
noise  powers  at  the  input  and  the  output  of  the  various  processors.  The  signal-only 
power  levels  are  estimated  by  substracting  the  noise-only  power  from  the  signal- 
plus-noise  levels.  When  the  2048  point  FFT  is  used,  the  average  noise-only  power  at 
198.73  Hz  is  71.979  dB  (126  snapshots)  at  the  input  of  the  processors,  while  the 
signal-plus-noise  power  is  84.054  dB  (126  snapshots).  Then,  the  input  SNR  is  on  the 
order  of  1 1.8  dB,  which  is  close  to  the  estimated  input  SNR  of  1 1  dB.  When  a  16 
point  FFT  is  used,  the  input  signal-to-noise  ratio  has  been  theoretically  reduced  by 
21  dB  and  should  be  on  the  order  of  -9.2  dB.  Because  of  the  very  short  FFT,  200  Hz 
signal  free  data  collected  18  minutes  after  the  signal  data  are  used  to  estimate 
noise  powers  at  the  input  and  output  of  the  beamformers.  The  average  noise  power 
at  the  input  of  the  processors  is  equal  to  93.6  dB  (126  snapshots)  while  the  signal- 
plus-noise  average  power  across  the  array  is  94.87  dB  (126  snapshots).  Using  these 
raw  data,  the  input  signal-to-noise  ratio  is  found  equal  to  -4.7  dB,  which  is  too  large 
compared  to  its  expected  value.  One  should  note  at  this  point  that  the  input  signal- 
to-noise  ratio  is  very  sensitive  to  small  errors  in  the  noise  and  signal  values.  For 
instance  an  increase  of  one  dB  in  noise  power  produces  an  input  SNR  of  -12  dB. 
Since  signal  and  noise  data  do  not  correspond  to  the  same  absolute  time, 
irreductible  errors  will  be  found  in  the  power  estimates.  The  array  gain  will  be 
estimated  in  the  following  with  signal-plus-noise  levels  from  data  starting  at 
22:02:49  GMT,  noise-only  levels  from  data  starting  at  22:20:17  GMT  and  an  input 
signal-to-noise  ratio  equal  to  -9.2  dB.  Table  4.2  summarizes  the  array  gain  values 
produced  by  each  method. 

The  values  of  the  array  gain  for  the  conventional  FFT  beamformer  should  be 
independent  of  the  input  signal-to-noise  ratio  and  the  1  dB  error  found  is  due  to 
the  erroneous  estimation  of  the  noise  powers  in  the  lower  input  SNR  case.  Given 
the  time  elapsed  between  the  signal-plus-noise  power  estimates  and  the  noise-only 
power  estimates,  such  an  error  can  be  considered  acceptable.  The  array  gains 
produced  by  full  aperture  and  subaperture  MV  processing  schemes  are  very  close  to 
the  conventional  FFT  beamformer  array  gain  in  the  case  of  the  -9.2  dB  input  SNR. 
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When  the  input  SNR  is  1 1.8  dB,  the  array  gain  varies  depending  on  the  method  used 
and  is  lower  than  the  conventional  FFT  beamformer  array  gain  by  at  least  a  few  dB. 
This  is  expected  because  of  mismatch.  For  both  input  SNRs,  the  experimental 
results  are  consistent  with  the  simulations  of  Chapter  3.  The  array  gains  obtained 
with  the  beamformers  on  real  or  simulated  data  have  the  same  order  and  have 
values  below  the  theoretical  array  gain  value  predicted  by  the  10  log  M  rule,  or  21 
dB. 


4.7.  Conclusions 

The  processing  schemes  proposed  in  Chapter  3  were  used  to  process  real  data. 
Resolution  improvements,  compared  to  conventional  FFT  beamforming,  are 
obtained  with  the  adaptive  MV  processor.  The  finer  structure  in  the  angular  spectra 
gives  some  insights  on  the  propagation  physics.  Using  the  MV  processor  on 
subaperture  indicates  that  the  pressure  field  is  highly  inhomogeneous  and 
multipath  arrivals  insonify  only  part  of  the  900  m  aperture.  These  complicated 
propagation  characteristics  were  succesfully  modelled  with  both  the  Generic  Sonar 
Model  and  the  ATLAS  normal  mode  model.  After  a  qualitative  study  of  the  various 
estimates  of  the  angular  spectra,  resolution  and  signal-plus-noise  to  noise  ratios  are 
quantified  for  a  high  and  medium  input  signal-to-noise  ratio  (11  dB  and  -9  dB).  The 
performance  indices  are  consistent  with  the  patterns  outlined  in  the  simulations  of 
Chapter  3. 


Array  Gain  (dB) 

Method 

Input  SNR 
11.8  dB 

Input  SNR 
-9.2  dB 

Bartlett 

14.9 

13.9 

MV  SSL  -  92 

8.0 

13.2 

MV  SSL  -  78 

12.8 

13.4 

MV  SSL  -  64 

12.9 

13.1 

MV  2  SSL  -  48 

8.4 

12.9 

MV  2  SSL  =  32 

11.2 

12.3 

Table  4.2:  Array  Gain. 
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Figure  4.1:  Sound  speed  profile  at  the  R/P  FLIP. 
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Figure  4,3:  Distribution  of  signal  and  noise  power  across  the  VLA. 
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Figure  4.4: 

200  Hz. 


Conventional  FFT  beamformer  output  for  signal  (left)  and  noise  (right)  at 
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Figure  4.5:  GSM  ray  traces. 
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Figure  4.6:  GSM  power  distribution  across  the  VLA. 


Figure  4.7:  GSM  eigenrays  at  162  km  and  at  the  sound  axis. 
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Figure  4.10:  Conventional  FFT  beamformer  output  on  the  ATLAS  field 


Figure  4.1 1:  Conventional  FFT  beamformer  output. 


Figure  4.12:  MV  processor  output  after  spatial  smoothing  on  64  sensors. 
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Figure  4.13:  MV  processor  output  after  modified  spatial  smoothing  on  64  sensors. 


Figure  4.14:  MV  processor  output  after  modified  spatial  smoothing  on  78  sensors. 


Figure  4.15:  MV  processor  output  after  modified  spatial  smoothing  on  92  sensors. 
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Figure  4.16:  Upper  subaperture  MV  spectrum  (modified  smoothing  length  of  48). 


Figure  4.17:  Lower  subaperture  MV  spectrum  (modified  smoothing  length  of  48). 


Figure  4.18:  Average  of  the  two  non-overlapping  subaperture  MV  angular  spectra. 
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Figure  4.19:  GSM  eigenrays  between  10°  and  15°. 
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Figure  4.20:  GSM  ray  traces  between  160  and  164  km. 
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Figure  4.21:  ATLAS  modal  eigenfunctions. 


Figure  4.22:  Upper  and  lower  subaperture  MV  spectra  for  a  single  mode  pressure  field. 


Figure  4.23:  Upper  and  lower  subaperture  MV  spectra  after  the  Snell’s  law  remapping. 
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Figure  4.24:  ATLAS  transmission  loss 
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25:  Upper  and  lower  subaperture  MV  spectra  of  the  ATLAS  field  at  165  km. 


Figure  4.26:  Average  of  the  subaperture  MV  angular  spectra  of  the  ATLAS  field. 
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resolution  performance  study:  1 1  dB  input  SNR. 


resolution  performance  study:  -9  dB  input  SNR. 


5.  Matched  Field  Processing  on  the  VLA  Data 


5.1.  Overview 

Based  on  experimental  data  collected  by  the  MPL  digital  array,  this  chapter 
demonstrates  that  matched  field  processing  (MFP)  is  a  viable  processing  technique 
at  200  Hz  for  a  900  m  aperture  array.  The  available  knowledge  of  the  environment 
facilitates  good  acoustic  modeling  which  leads  to  the  detection  and  localization  of  a 
source  at  20  m  depth  and  165  km  range.  The  Minimum  Variance  processor 
performs  poorly  as  a  result  of  mismatch  due  to  a  1°  array  tilt  (which  is  severe  at 
200  Hz).  The  array  gain  estimates  obtained  with  the  MFP  Bartlett  processor  at  high 
(11.8  dB)  and  low  (-9.2  dB)  input  signal-to-noise  ratio  are  higher  than  the  array  gain 
estimates  of  the  conventional  FFT  beamformer.  Based  on  these  experimental  results, 
matched  field  processing  is  shown  to  be  a  better  beamforming  approach  for  signal 
detection  than  is  conventional  processing. 


5.2.  Introduction  to  Matched  Field  Processing 

A  conventional  beamformer  which  processes  the  data  from  a  vertical  line  array 
to  estimate  an  angular  spectrum,  addresses  both  the  detection  and  estimation 
problems.  The  first  function  of  a  beamformer  is  to  detect  sound  sources  such  as 
ships  by  taking  advantage  of  the  spatial  processing  gain.  The  detection  problem  is 
solved  by  inspecting  the  peaks  of  the  angular  spectrum.  These  same  peaks  also 
allow  the  estimation  of  the  directions  of  arrival.  A  vertical  line  array  processed  in 
the  conventional  fashion  does  not  provide  a  direct  answer  to  the  localization 
problem.  To  localize  the  sound  sources,  one  has  to  compare  the  experimental 
vertical  angular  spectrum  with  a  predicted  arrival  structure  which  assumes  a  known 
acoustic  environment  and  source  position.  This  "matched  filtering"  approach,  which 
is  implemented  indirectly  through  the  use  of  conventional  beamformers,  is  carried 
out  directly  by  matched  field  processing  (MFP)  [Bucker;  1976,  Baggeroer;  1988, 

Fizell;  1987,  Smith;  1989]. 

Instead  of  assuming  a  simple  propagation  model  involving  planar  and 
uncorrelated  arrivals,  matched  field  processing  is  based  on  a  more  general  model 
which  includes  wavefront  curvature  and  multipath  arrivals.  It  performs  a  global 
search  to  detect  and  localize  signals  by  relating  a  single  pressure  contribution  to 
each  source.  The  model  array  covariance  matrix  for  k  uncorrelated  sources  is 

R=  sVd.d^  +  Q  (5.1) 

I  m  0 

where  Q.is  the  noise  covariance  matrix,  at2  is  the  power  of  the  ith  signal.  The 
normalized  vector  df  (d,H  d,  =  1)  does  not  contain  simple  plane  wave  phase 
relationships  indexed  on  a  local  angle  as  in  conventional  beamforming,  but  is  the 
total  pressure  field  corresponding  to  the  ith  source.  Matched  field  processing  uses  a 
synthetic  pressure  field  corresponding  to  a  particular  source  range  and  depth  as  a 
replica  vector.  Direct  localization  is  obtained  from  the  resulting  ambiguity  surface 
peaks,  since  the  search  space  is  indexed  on  source  range  and  depth.  Matched  field 
processing  involves  complex  modeling  in  order  to  compute  the  array  manifold  but 
uses  the  same  processing  structures  as  conventional  beamforming  to  measure  the 
correlation  between  the  data  and  the  predicted  pressure  fields.  For  this  reason, 
matched  field  processors  are  called  generalized  beamformers. 

The  MFP  signal  model,  which  includes  wavefront  curvature  and  multipath 
arrivals,  is  closer  to  reality  than  the  simple  signal  model  of  the  conventional 
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beamformers  discussed  in  the  previous  chapters.  This  suggests  that  matched  field 
processing  produces  a  larger  array  gain  and  is  a  better  detection  tool  than  the 
conventional  beamformers,  in  addition  to  providing  direct  source  localization.  The 
array  gain  achieved  by  coherently  processing  M  sensors  is  10  log  M  [Pillai;  1989] 
assuming  coherent  signals  and  incoherent  noise.  Such  array  gain  values  cannot  be 
expected  by  the  conventional  beamformers  since  multipath  arrivals  are  searched  in 
angle,  as  if  they  were  uncorrelated.  Matched  field  processing,  which  uses  the  full 
wavefield  and  recombines  the  energy  of  all  the  multipaths  to  calculate  the 
correlation  between  the  measured  and  the  predicted  field,  theoretically  will  yield 
this  10  log  M  array  gain  [Heitmeyer;  1985,  Yang;  1989].  The  issue  of  the  array  gain 
achieved  by  conventional  beamformers  such  as  the  FFT  beamformer,  was  addressed 
in  simulations  in  Chapter  3  and  on  real  data  in  Chapter  4.  The  simulations 
presented  in  Chapter  3  involved  a  100  Hz  CW  source  deployed  at  5  m  depth  and  at 
56  km  range  from  a  128  sensor  VLA  cut  for  107  Hz.  A  realistic  ocean  waveguide  was 
considered  with  a  sound  speed  profile  characteristic  of  the  North-East  Pacific.  The 
matched  field  signal  model  was  used  to  produce  simulated  data  using  a  normal 
mode  model.  The  synthetic  wavefield,  then,  was  beamformed  with  the  conventional 
beamformers  (the  FFT  beamformer  and  the  Minimum  Variance  beamformer).  The 
array  gain  was  computed  and  found  to  at  most  14  dB,  i.e.  7  dB  below  the  expected 
array  gain  of  10  log  M  or  21  dB  for  M  =  128.  Array  gain  estimates  on  real  data  in 
Chapter  4  were  found  to  be  at  most  14  to  15  dB  on  200  Hz  CW  signals  recorded  by 
the  MPL  digital  acoustic  array,  a  120  sensor  vertical  line  array,  also  with  a  21  dB 
expected  array  gain. 

The  ability  to  detect  and  localize  with  the  matched  field  processing  techniques 
is  highly  dependent  on  the  knowledge  of  the  environment  and  the  subsequent 
acoustic  modeling  involved  in  the  computation  of  the  array  manifold  [Porter;  1987, 
Del  Bazo;  1988,  Feuillade;  1989,  Tolstoy;  1989,  Baxley;  1989].  Matched  field 
processing  requires  a  detailed  knowledge  of  the  sound  speed  structure  of  the 
waveguide,  its  eventual  range  dependence  as  well  as  other  environmental 
parameters.  An  erroneous  modeling  produces  mismatch  which  can  severly  degrade 
the  matched  field  processing  performances.  In  comparison,  conventional 
beamformers  require  no  or  little  knowledge  of  the  environment.  The  curved 
wavefront  beamformers.  proposed  in  Chapter  2  and  3,  only  requires  knowledge  of 
the  sound  speed  across  the  array.  An  issue  with  matched  field  processing  is  to 
determine  whether  or  not  a  sufficient  knowledge  of  the  environment  is  available 
before  attempting  to  process  real  data. 

Good  environmental  information  was  collected  during  the  September  1987  VLA 
Experiment  and  allowed  an  accurate  acoustic  modeling  of  the  data  of  Chapter  4. 
Matched  field  processing  is  attempted  here  on  the  same  data  set.  Localization  and 
processing  gain  performances  are  examined. 


5.3.  Matched  Field  Processing  on  Real  Data 


5.3.1.  Description  of  the  Processing 

Matched  field  processing  techniques  blend  signal  processing  and  acoustic 
modeling  as  shown  in  Figure  5.1.  These  two  areas  were  investigated  in  Chapter  4. 
After  data  quality  checks,  the  array  covariance  matrix  was  generated  and  processed 
by  the  beamformers  to  produce  angular  spectra.  Acoustic  modeling  with  the  ATLAS 
normal  mode  model  has  been  performed  using  the  available  environmental  data, 
compared  to  the  experimental  results  and  is  now  used  to  create  the  MFP  array 
manifold.  The  matched  field  techniques,  then,  process  the  data  covariance  matrix 
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with  the  set  of  replica  vectors  generated  by  the  normal  mode  model,  using  the 
conventional  or  Bartlett  processor  and  the  MV  processor. 

The  same  data  segments,  as  in  the  case  of  the  conventional  beamformers  of 
Chapter  4,  are  used  to  estimate  the  data  covariance  matrix.  The  covariance  matrix  is 
obtained  by  averaging  126  snapshots  computed  with  2048  point  FFTs  to  get  a  11.8 
dB  input  SNR  and  or  a  16  point  FFT  to  obtain  a  -9.2  dB  input  SNR.  Both  signal-plus- 
noise  and  noise-only  data  are  processed  by  the  MFP.  When  a  2048  point  FFT  is  used, 
the  signal-plus-noise  and  noise-only  data  correspond  to  the  same  absolute  time  but 
have  slightly  different  frequencies.  The  signal  FFT  bin  is  #  1844  or  200.2  Hz,  and 
the  noise  FFT  bin  is  five  FFT  bins  away  (#  1839)  or  198.73  Hz.  In  the  case  of  the  16 
point  FFT,  the  bin  width  is  on  the  order  of  30  Hz  and  one  uses  noise-only  data 
recorded  18  minutes  after  the  transmissions,  at  the  same  center  frequency  (Section 
4.6).  The  array  covariance  matrix,  for  both  input  SNR,  can  be  easily  inverted  with  a 
small  stabilization  factor  as  in  Chapter  4,  but  will  be  stabilized  here  with  a  fraction 
of  noise  7  equal  to  0.001,  0.01  and  0.1  to  enhance  the  robustness  of  the  MV 
processor  (Section  1.6.3. 2). 

The  repiica  vectors  are  computed  for  a  source  in  a  spatial  window  extending  in 
range  from  152.5  to  177.5  km,  and  in  depth  from  5  and  600  m.  The  sampling  rates  in 
range  and  in  depth  are  respectively  250  m  and  5  m.  The  sound  speed  profile  is 
derived  from  the  CTD  and  XBT  data  recorded  at  the  R/P  FLIP  at  the  time  of  the 
experiment.  The  bottom  is  characterized  by  a  GSM  Bottom  Province  Type  3,  its 
influence  is  weak  because  the  propagation  experiment  takes  place  in  deep  water. 
The  ambiguity  surfaces  computed  by  the  Bartlett  and  MV  processors  are  plotted  in 
logarithmic  gray  levels.  In  each  case,  the  gray  levels  span  the  full  dynamic  range  of 
the  surface  which  is  normalized  to  its  highest  peak.  The  peak  value  is  recorded  and 
normalized  to  yield  power  in  dB  re  /xPa. 


5.3.2.  11.8  dB  Input  SNR  Case 

The  signal  with  11.8  dB  input  SNR  is  processed  by  the  matched  field  processor 
and  the  results  presented  in  Figure  5.2  and  5.3.  The  source  at  20  m  is  detected  at 
the  expected  source  location  as  shown  in  Table  5.1. 


Processor 

Highest  Peak 
Power 
dB 

Highes 

Loca 

Range  (km) 

t  Peak 
tion 

Depth  (m) 

Bartlett 

79.0 

163.50 

25 

49.6 

164.25 

20 

MV  (7  =  0.01) 

54.4 

164.25 

20 

MV  (7  =  0.1) 

60.0 

164.25 

20 

Table  5.1:  Matched  field  processing  (11.8  dB  input  SNR). 


The  ambiguity  surfaces  are  granular  with  a  large  number  of  sidelobes.  One 
recognizes  in  the  Bartlett  ambiguity  surface  the  usual  pattern  which  is  similar  to  a 
field  plot  at  a  convergence  zone.  Using  the  MV  reduces  these  sidelobes,  which  have 
less  and  less  structure  as  the  fraction  of  white  noise  7  increases.  The  dynamic  range 
in  the  ambiguity  surfaces  of  the  MV  processor  is  small  compared  to  the  one  of  the 
Bartlett  processor.  The  MV  processor  suffers  large  loss  due  to  mismatch  in  this 
high  signal-to-noise  ratio  case.  As  the  amount  of  stabilization  in  the  covariance 
matrix  increases,  the  robusness  is  increased  and  the  highest  peak  absolute  level 
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increases. 

Matched  field  processing  on  the  noise-only  data  are  presented  in  Figures  5.4 
and  5.5.  The  ambiguity  surfaces  indicate  darker  areas  below  300  m,  deep  in  the 
water  column.  The  highest  peak  appears  at  large  depths  at  440  m  or  at  550  m.  The 
MV  processor  also  indicates  a  surface  source  at  the  source  range.  Since  a  Kaiser- 
Bessel  window  with  57  dB  sidelobe  rejection  level  is  used  in  the  time  to  frequency 
transform,  this  is  not  the  results  of  200  Hz  signal  leakage.  This  peak  in  the 
ambiguity  surface  probably  correspond  to  the  radiated  noise  of  the  support  ship 
from  which  the  source  was  deployed.  It  could  also  be  a  sidelobe  from  a  distant 
(unknown)  source.  The  noise  data  are  processed  with  198.7  Hz  replica  vectors, 
although  using  the  200  Hz  replica  vectors  of  the  signal-plus-noise  case  would 
produce  negligible  differences  in  the  ambiguity  surface  structure  and  peak  levels. 
On  the  other  hand,  simulations  involving  a  source,  in  a  spatially  white  noise  with  10 
dB  input  SNR  show  that  the  mismatch  due  to  a  frequency  offset  of  2  Hz  has  a  major 
impact  on  the  MV  processor  with  mismatch  losses  on  the  order  of  tens  of  dB. 


5.3.3.  -9.2  dB  Input  SNR  Case 

The  results  for  the  low  input  SNR  case  are  presented  in  Figures  5.6  and  5.7  for 
the  signal  data  and  in  Figures  5.8  and  5.9  for  the  noise-only  data.  The  source  is 
detected  at  its  correct  location  by  the  Bartlett  processor.  The  conventional 
processor  ambiguity  surface  for  signal  data  is  similar  to  the  high  input  SNR  case. 
Some  structure  appears  in  the  surface  below  a  depth  of  300  m  and  corresponds  to 
the  noise  field.  Below  300  m.  a  great  similarity  is  observed  between  the  signal-plus- 
noise  and  noise-only  ambiguity  surfaces.  The  highest  peak  in  the  ambiguity  surface 
produced  by  the  Minimum  Variance  processor  does  not  corresponds  to  the 
expected  source  (20  m  depth  and  165.5  km  range).  Nevertheless,  a  peak  can  be 
found  next  to  the  expected  source  position  at  a  depth  of  20  m  and  a  range  of 
166.25  km,  with  a  power  level  close  to  the  highest  power  level.  This  signal  becomes 
better  defined  as  the  amount  of  stabilization  in  the  array  covariance  matrix  is 
increased  and  will  be  used  in  the  performance  study.  Table  5.2  summarizes  the 
signal  information  with  the  highest  peak  power  level  and  location  as  well  as  the 
probable  source  strength  and  location,  based  on  the  ambiguity  surfaces.  Matched 
field  processing  on  the  noise-only  data  which  were  recorded  18  minutes  after  the 
signal  data  produces  ambiguity  surfaces  which  have  a  similar  noise  structure  as  in 
the  high  signal-to-noise  ratio.  The  ambiguity  surfaces  have  darker  areas  below  300 
m  in  the  deep  sound  channel.  Since  the  input  SNR  is  lower  by  20  dB  compared  to 
the  previous  section,  no  signal  is  detected  at  the  support  ship  range. 


Processor 

H 

ighest  Peak 

Probable  Target 

Level 

(dB) 

Range 

(km) 

Depth 

(m) 

Level 

(dB) 

Range 

(km) 

Depth 

(m) 

Bartlett 

80.97 

166.5 

20 

80.97 

166.5 

20 

MV  h  =  0.001) 

62.85 

156.25 

220 

61.54 

166.25 

20 

MV  (7  =  0.01) 

65.56 

161.75 

190 

65.33 

166.25 

20 

o 

II 

> 

2 

70.64 

165.25 

95 

70.61 

166.25 

20 

Table  5.2:  Matched  field  processing  (-9.2  dB  input  SNR). 
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5.4.  Matched  Field  Processing  Simulations 

After  these  qualitative  experimental  results,  which  show  excellent  localization 
performance,  simulations  are  now  presented  in  an  effort  to  understand  the 
experimental  ambiguity  surfaces. 


5.4.1.  No  Mismatch  Simulations 

A  perfect  situation  is  now  assumed  where  there  is  no  noise  and  no  mismatch 
(i.e.  the  replica  vectors  computed  by  the  ATLAS  model  constitute  the  "truth").  A  200 
Hz  source  is  simulated  at  a  range  of  165  km  and  a  depth  of  20  m.  The  ambiguity 
surfaces  for  the  Bartlett  processor  and  the  three  different  implementations  of  the 
MV  processor  are  computed  and  plotted  in  Figures  5.10  and  5.1 1.  The  source  is 
correctly  localized.  Since  there  is  no  mismatch  and  no  noise,  the  dynamic  range  of 
the  ambiguity  surfaces  is  much  larger  than  in  the  case  of  the  real  data.  The  Bartlett 
processor  produces,  in  all  these  simulation  cases,  a  large  number  of  sidelobes.  It 
exhibits  the  hat  shape,  also  observed  in  the  real  data  ambiguity  surfaces.  The  MV 
processor  has  a  great  sensitivity  to  mismatch  in  this  noise-free  simulations  and 
produces  a  single  peak  in  the  range-depth  cell  of  the  simulated  source.  These 
simulation  results  indicate  that  there  is  some  mismatch  between  the  replica  vectors 
and  the  data. 


5.4.2.  Analysis  of  Mismatch 

There  are  several  potential  causes  for  the  large  mismatch  observed  in  the 
matched  field  processing  results.  Although  sound  speed  information  at  the  R/P 
FLIP  is  very  good,  a  range  independent  medium  is  assumed  which  may  not  be 
accurate  for  ranges  on  the  order  of  165  km.  There  is  no  available  range  dependent 
information  to  improve  the  processing.  As  mentioned  earlier,  the  bottom  should 
have  only  a  weak  influence  in  the  processing  results  since  this  is  a  deep  water 
environment.  Another  cause  for  mismatch  is  the  position  errors  of  the  array 
sensors.  The  VLA  array  was  deployed  in  a  three  transponder  network  and  one 
channel  per  array  section  was  navigated  during  the  experiment  [Sotirin;  1989]. 
Unfortunately  navigation  data  are  not  available  during  the  data  segment  processed 
here  (Tape  918)  because  of  temporary  hardware  problems.  Navigation  information 
is  available  for  Tape  915,  40  minutes  earlier  and  is  summarized  in  Figure  5.12.  The 
array  shapes  are  plotted  in  the  North-South  and  East-West  vertical  planes.  These 
typical  navigation  data  indicate  that  the  array  behaves,  to  the  first  order,  like  a 
pendulum  and  the  largest  horizontal  displacement  across  the  array  aperture  is 
bounded  by  15  m.  Such  displacement,  corresponding  to  a  1“  tilt,  represents  two 
wavelengths  at  200  Hz  and  is  likely  to  have  a  major  impact  on  the  matched  field 
processing  results.  This  source  of  mismatch  is  now  investigated. 


5.4.3.  Tilted  Array  Simulation 

A  tilted  line  array  is  simulated  with  horizontal  displacement  of  15  m.  The 
source  is  assumed  at  165  km  from  the  array  top  element  and  at  a  depth  of  20  m. 
This  tilted  array  matched  field  processing  simulation  is  summarized  in  Figure  5.13. 
The  replica  vectors,  which  were  used  earlier,  assume  a  straight  line  vertical  array. 
The  impact  of  mismatch  due  to  array  tilt  is  investigated  in  a  noise  free  case,  and  the 
ambiguity  surfaces  are  plotted  in  Figure  5.14  and  5.15.  The  matched  field 
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processing  results  also  are  summarized  in  Table  5.3.  Array  tilt  results  in  enough 
mismatch  to  significantly  reduce  the  maximum  peak  power  and  shift  the  position  of 
the  peak  in  range  and  depth.  The  mismatch  is  sufficient  to  reduce  in  a  dramatic  way 
the  signal-to-background-noise  ratio.  The  estimated  position  of  the  source  is  not 
the  true  one  and  is  offset  by  4  km  in  range.  These  simulated  results  suggest  that 
array  tilt  can  be  held  responsible  for  the  large  mismatch  observed  in  the  real  data 
ambiguity  surfaces. 


Bartlett 

Proce: 

MV 

7  =  lO"3 

;sor 

MV 

7  =  lO'2 

MV 

7  =  lO"1 

Range  of  Max. 

160.75  km 

161  km 

161  km 

161  km 

Surface  Max. 

-8.7  dB 

-49.8  dB 

-39.8  dB 

-29.8  dB 

Table  5.3:  Tilted  array  matched  field  processing  simulation 


5.5.  Matched  Field  Performance  Study 

In  addition  to  the  localization  performance  criteria,  one  is  interested  in 
performance  measures  that  characterize  matched  field  processing.  Several 
performance  criteria  have  been  proposed  in  the  literature,  such  as  the  relative  peak 

P  _  n 

to  mean  background  level  given  by - ^  where  P  is  the  peak  level  and  n  the  mean 

P 

level  [Del  Bazo;  1989],  or  the  peak  level  above  the  mean  level  relative  to  the 
standard  deviation  of  the  ambiguity  surface  <r, - =-  [Feuillade;  1989],  or  the  peak 

(T 

relative  to  the  highest  sidelobe  level.  None  of  these  common  performance  criteria  is 
satisfying  because  they  do  not  fully  characterize  the  ambiguity  surface  and,  given 
the  complexity  of  the  experimental  ambiguity  surfaces,  they  would  indicate  poor 
quality  results  most  of  the  time.  In  the  following,  ihe  emphasis  is  put  on  the  array 
gain  or  processing  gain  produced  by  matched  field  processing.  A  main  advantage 
of  using  array  gain  as  a  performance  index  is  that  it  allows  an  intercomparison 
between  matched  field  processing  and  the  conventional  beamformers.  The  array 
gain  is  estimated  as  in  Chapter  4.  The  peak  indicating  the  signal  (generally  the 
highest  level  peak)  is  used  to  measure  the  output  signal-plus-noise  level.  The 
output  noise  level  is  measured  for  the  same  range  and  depth  cell  from  the  noise- 
only  ambiguity  surfaces.  The  array  gain  can  then  be  computed  using  the  estimate 
of  the  input  signal-to-noise  ratio. 


5.5.1.  11.8  dB  Input  SNR  Case 

Signal-plus-noise  and  noise-only  data  are  summarized  in  Table  5.4  together 
with  the  array  gain.  The  array  gain  of  the  Bartlett  processor  is  on  the  order  of  16  dB. 
1  dB  above  the  array  gain  achieved  by  the  conventional  FFT  beamformer  in  Chapter 
4,  on  the  same  data  set.  This  computed  array  gain  is  lower  than  expected  since  it  is 
5  dB  below  the  21  dB  maximum  expected  array  gain.  The  MV  processor  is  very 
sensitive  to  mismatch  at  high  signal-to-noise  ratios  and  its  processing  gain  is  poor. 

It  increases  with  the  covariance  matrix  stabilization  thus  indicating  serious 
problems. 
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Processor 

S+N 

dB 

N 

dB 

AG 

dB 

Bartlett 

79 

51.47 

15.7 

MV  (7  =  0.001) 

49.6 

36.1 

1.5 

MV  (7  =  0.01) 

54.4 

39.9 

2.5 

MV  (7  =  0.1) 

60 

45.3 

2.7 

Table  5.4:  Array  gain  calculation  (1 1.8  dB  input  SNR) 


5.5.2.  -9.2  dB  Input  SNR  Case 

The  probable  target  is  used  to  measure  the  signal-plus-noise  peak  levels.  The 
MFP  output  for  noise-only  data  are  computed  at  the  probable  source  range  and 
depth  for  all  four  processors.  The  signal-plus-noise  and  noise-only  levels  are 
summarized  in  Table  5.5.  The  signal-only  power  levels  are  estimated  and  the  output 
signal-to-noise  ratios  and  array  gains  computed.  The  array  gain  for  the  Bartlett 
processor  is  on  the  order  of  17.5  dB  which  is  higher  than  the  array  gain  for  the  high 
input  SNR  case.  The  presence  of  a  peak  in  the  ambiguity  surface  at  the  source 
location  even  with  the  source  turned  off  explains  perhaps  such  a  discrepency.  The 
MFP  Bartlett  processor  has  an  array  gain  3  dB  higher  than  the  conventional  FFT 
beamformer.  Because  of  the  lower  input  SNR,  the  MV  processor  is  less  sensitive  to 
mismatch  and  performs  more  reasonably  than  in  the  high  input  SNR  case.  The  MV 
array  gain  is  between  11  and  13  dB  and  decreases  as  the  amount  of  stabilization  in 
the  covariance  matrix  increases.  This  is  consistent  with  the  fact  that  additional 
artificial  noise  is  added  in  the  process. 


Processor 

S+N 

dB 

N 

dB 

AG 

dB 

Bartlett 

80.97 

72.13 

17.43 

61.54 

56.23 

13.0 

Efi92B§X§Hl 

65.33 

60.86 

12.35 

70.61 

66.58 

11.05 

Table  5.5:  Array  gain  calculations  (-9.2  dB  input  SNR) 


5.6.  Conclusions 

This  chapter  shows  that  matched  field  processing  is  a  viable  processing 
technique  at  200  Hz  with  a  900  m  long  aperture  array.  The  source  at  20  m  depth 
and  165  km  range  is  detected  at  two  input  signal-to-noise  ratios  equal  to  1 1.8  and 
-9.2  dB.  The  source  localization  is  excellent  in  both  range  and  depth.  The 
environmental  information  collected  at  the  time  of  the  experiment  at  the  R/P  FLIP 
allows  an  appropriate  modeling  of  the  acoustic  transmissions.  Mismatch  exists 
because  of  a  severe  1°  array  tilt.  The  lack  of  robustness  of  the  MV  processor  is 
observed  with  these  real  data.  The  estimated  array  gain  of  matched  field  processing 
using  the  conventional  processor  is  found  larger  then  the  processing  gain  estimates 
of  the  conventional  FFT  beamformer.  Even  under  large  mismatch  conditions, 
matched  field  processing  appears  to  be  a  better  beamforming  approach  for  signal 
detection  than  is  conventional  processing. 
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Figure  5.1:  Matched  Field  Processing. 
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Figure  5.3:  Matched  field  processing  results  on  signal,  pan  2,  (11.8  dB  input  SNR  case). 
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Figure  5.5:  Matched  field  processing  results  on  noise,  pan  2,  (1 1.8  dB  input  SNR  case). 
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Figure  5.8:  Matched  field  processing  results  on  noise,  part  1,  (-9.2  dB  input  SNR  case). 
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Figure  5.10:  Matched  field  processing  simulation,  part  1,  (no  mismatch,  no  noise). 
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Figure  5.12:  Array  navigation  data  for  Tape  915. 
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Figure  5.13:  Tilted  array  simulation  geometry. 
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Figure  5.14:  Matched  field  processing  simulation,  part  1,  (tilted  array). 
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Figure  S.15:  Matched  field  processing  simulation,  part  2,  (tilted  array). 


Conclusions 


Various  approaches  to  process  the  data  from  large  aperture  vertical  line 
arrays  have  been  examined  in  this  dissertation.  Depending  on  the  available 
knowledge  of  the  environment,  beamforming  or  matched  field  processing  can  be 
used  to  enhance  the  detection  of  signals  through  an  increase  in  signal-to-noise 
ratio,  and  also  provide  direction  of  arrival  estimation  or  range-depth  localization. 

In  a  scenario  where  little  is  known  on  the  oceanic  environment,  beamform¬ 
ing  is  performed  to  detect  signals  and  estimate  the  field  vertical  arrival  structure. 
The  use  of  the  adaptive  Minimum  Variance  processor,  with  large  aperture  vertical 
line  arrays,  has  been  studied  and  its  performance  compared  to  the  conventional  FFT 
beamformer.  As  a  preliminary  to  processing  real  data,  the  signal  cancellation  prob¬ 
lem  of  correlated  arrivals  is  examined  and  the  spatial  smoothing  techniques  inter¬ 
preted.  The  mismatch  due  to  wavefront  curvature  also  is  evaluated.  It  is  shown  that 
wavefront  curvature  should  be  taken  into  account  for  a  3000  m  aperture  array  and 
for  arrivals  with  low  angle  with  respect  to  the  horizontal.  On  the  other  hand,  plane 
wavefront  replica  vectors  can  be  used  with  a  900  m  aperture  array  even  at  high 
input  signal-to-noise  ratio  (on  the  order  of  10  dB),  when  there  is  no  arrival  shal¬ 
lower  than  ±10°.  The  estimation  of  a  100-sensor  array  covariance  matrix  requires  a 
large  number  of  data  snapshots  from  a  potentially  nonstationary  pressure  field.  For 
this  reason,  subaperture  Minimum  Variance  processing  was  proposed.  Using 
subapertures  reduces  the  computational  load  in  addition  to  the  number  of 
snapshots  entering  into  the  covariance  matrix,  but  at  the  cost  of  processing  gain 
and  resolution  because  of  the  reduction  of  the  aperture  length.  The  full  and 
subaperture  processing  schemes  are  validated  in  simulations  and  used  on  real  data 
collected  by  the  MPL  digital  array.  The  multipath  arrivals  are  found  to  illuminate 
only  parts  of  the  array  and  based  on  the  particular  data  set  presented  in  this  disser¬ 
tation,  a  real  wavefield  may  be  highly  inhomogeneous,  Subaperture  MV  processing 
is  shown  to  be  valuable  in  illustrating  the  great  multipath  variability,  but  averaging 
subaperture  angular  spectra  may  not  enhance  the  peaks  in  the  composite  estimate. 
The  processing  gain  obtained  with  the  MV  processor  is  lower  than  the  conventional 
FFT  beamformer  because  of  spatial  smoothing  or  residual  signal  cancellation,  and 
mismatch.  On  the  other  hand,  resolution  can  be  significantly  improved,  especially 
for  high  level  arrivals. 

The  array  gain  estimates  of  the  conventional  FFT  beamformer  were  found 
between  14  and  15  dB  in  the  simulation  and  real  data  cases,  i.e.  6  to  7  dB  lower  than 
what  matched  field  processing  would  achieve  in  the  case  of  no  mismatch  (21  dB  or 
10  log  M  where  M  is  the  number  of  sensors).  Matched  field  processing  was 
attempted  on  200  Hz  real  data  signals  that  propagated  to  the  MPL  array  over  165 
km.  The  available  knowledge  of  the  environment,  the  acoustic  modeling  and  the 
good  quality  of  the  data  make  this  attempt  successful.  Generally,  the  source  is 
correctly  localized  in  range  and  depth.  The  main  cause  of  mismatch  is  identified  to 
be  the  array  tilt,  at  most  equal  to  1°  (which  corresponds  to  a  horizontal  displace¬ 
ment  of  two  wavelengths  at  200  Hz  from  the  top  to  the  bottom  of  the  array). 
Although  mismatch  is  present,  the  processing  gain  estimates  of  the  MFP  Bartlett 
processor  are  1  to  3  dB  larger  than  the  conventional  FFT  beamformer  array  gain 
estimates.  The  MV  processing  is  found  to  perform  poorly  as  a  result  of  mismatch 
and  lack  of  robustness. 

A  few  aspects  in  this  dissertation  should  be  the  subject  of  further  investiga¬ 
tions.  The  beamforming  schemes  proposed  here  can  be  implemented  with  proces¬ 
sors  other  than  the  Minimum  Variance  processor  which  is  known  to  be  very  sensi¬ 
tive  to  mismatch  at  high  signal-to-noise  ratio.  The  adaptive  procesor  under  addi¬ 
tional  constraints  (for  example  the  multiple  point  contraints  or  the  derivative 


constraints  [Van  Veen  1989]),  where  "real  world”  perturbations  a  accommodated, 
should  enjoy  more  robustness  which  is  an  important  characteristic  when  processing 
real  data.  The  problem  of  the  variance  of  the  angular  spectra  produced  by  the 
Minimum  Variance  processor  was  not  addressed  in  this  dissertation  and  should  be 
studied  in  order  to  treat  issues  related  to  detectability.  In  this  dissertation,  the 
emphasis  was  put  on  a  900  m  aperture  and  the  MPL  digital  array.  With  larger  aper¬ 
ture  arrays  expected  in  the  future,  it  would  be  interesting  to  test  further  via  simula¬ 
tion  and  the  processing  of  real  data  the  curved  wavefront  beamforming  techniques. 
Finally,  although  matched  field  processing  did  perform  reasonably  well,  efforts 
should  be  directed  toward  using  more  robust  processors  than  the  Minimum  Vari¬ 
ance  processor.  In  addition,  the  array  navigation  data  (sensor  position)  should  be 
incorporated  into  the  processing  and  the  increase  in  array  gain  measured. 
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